BBC  F1LE_CO0Q  ADA092918 


A  GENERAL  SURFACE  REPRESENTATION 
MODULE  DESIGNED  FOR  GEODESY 


.  HANS  -.SUNK-EL 


THE  OHIO  STATE  UNIVERSITY 
RESEARCH  FOUNDATION 


1 1 

JUNE  1980 

SCIENTIFIC  REPORT  NO.  3 

APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


AIR  FORCE  GEOPHYSICS  LABORATORY 
AIR  FORCE  SYSTEMS  COMMAND 
UNITED  STATES  AIR  FORCE 
HANSCOM  AFB,  MASSACHUSETTS  01731 


Me- 


1 1  005 


Unclassif  ed 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  ( When  Data  Entered) 


REPpRI  DOCUMENTATION  PAGE 

I^JWEP£LRX-N«^-a.E  Rl  I  9J - h.  GOV 

na  \  arn t  ifrvr^Q  d  diri/t  *'/  a  n 


70)  AFGL/lTR-8^-/Z)3;4 


3  A  r  P  READ  INSTRUCTIONS 

_ BEFORE  COMPLETING  FORM 

2.  GOVT  ACCESSION  NO.  3.  RECIPIENT'S  CATALOG  NUMBER 


I  A.  TITLE  (and  Subtitle) 


5.  TYPE  OF  REPORT  &  PERIOD  COVERED 


0  Y&H& NERAL  SURFACE  REPRESENTATION  ')  j 

MODULE  DESIGNED  FOR  GEODESY.  ^  o^r 

j  ^  ^  6‘  RP|po0rt '^o°.RG'2?£^T  NUMBER 


7.  AU^HgaW-'T 


Hans/Stfrikel 


li)^l 


|  9k  performing  organization  name  ano  aodress  / 

Department  of  Geodetic  Science  ,  \ 

The  Ohio  State  University  ^  V 

r.nlnmKn  c  f)Kin  4.^7  10  » 


8.  £°!jTR*CT  OQ-gRANT  NUMBERS 

)  F  196 28- 79 -C  -^75  ^ yS 

10.  PROGRAM  ELEMENT.  PROJECT.' TASK 
AREA  4  WORK  UNIT  NUMBERS 


II.  CONTROLLING  OFFICE  NAME  ANO  AOORESS  12  aFPfiaT  ntTF 

Air  Force  Geophysics  Laboratory  /7J)  J une.498^ 

Hanscom  AFB,  Massachusetts  017  31  (^/Uar-wuMULH’UF -pages  - 
_ Contract  Monitor  -  Bela  Szabo/LW _ 16  3  Pages _ 

IT  MONITORING  AGENCY  NAME  4  AQQR£SSaii_dlla«fiuil~i/am-Ca/maJlml-0»lci>)  IS.  SECURITY  CLASS,  (ol  this  report, 

lli>\  A  /  /  A  oA  /  v-r*rCA/-r-r  <<  )  Unclassified 


IfyMS-AtLSttZtf TJ  txt-i 


I  16.  DISTRIBUTION  STATEMENT  (ol  title  Report) 


is*,  declassification,  downgrading 
SCHEDULE 


Approved  for  public  release;  distribution  unlimited 


17.  DISTRIBUTION  STATEMENT  (ol  the  abstract  entered  In  Block  20,  II  dllletent  from  Report ) 


is.  supplementary  notes 


Tech,  Other 


*9  KEY  WORDS  (Continue  on  reverse  side  if  necessary  and  identify  by  block  number) 


Spline  representation,  Graphical  representation,  Contouring, 
Prediction,  Collocation 


20  ABSTRACT  i Continue  on  reverse  side  if  necessary  and  Identify  by  block  number) 

"  ^  Technological  developments  during  the  last  decade  have  provided 
geodesists  with  the  capability  of  obtaining  enormous  amounts  of 
data  on  and  outside  the  surface  of  the  earth.  In  order  to  understand 
the  general  behavior  of  the  data,  a  graphical  representation  of  the 
data  is  useful.  This  type  of  representation  can  be  obtained  by  using 
GSPP,  Geodetic  Science  Plotting  Package,  a  set  of  FORTRAN  rV 
_ subroutines.  This  report  describes  the  various  algorithms  of  GSPP 

DD  i  jan  73  1473  Unclassified  yt Lj  ^  ✓ 


Unclassified 


SECURITY  CLASSIFICATION  OF  THIS  PAGE(IVh»n  Data  Entarad) 


used  to  predict,  by  a  number  of  different  methods,  a  bicubic  spline 
'  surface  from  irregularly  distributed  data  and  to  produce  profiles, 
contour  maps  ard/or  3-*D  views  of  the  surface,  its  first  or  second 
derivatives.  In  addition,  complete  labelling  capabilities  are  included. 
Unique  features  of  GSPP  are  that  regions  in  which  contours  are  to  be 
included /excluded  can  be  defined  and  that  contour  maps  can  be  pro¬ 
duced  according  to  a  user  defined  projection. 


Unclassified 


SECURITY  CLASSIFICATION  OF  THIS  PACEfHTun  Data  Entarad) 


iii 


FOREWORD 

This  report  was  prepared  by  Dr.  Hans  Siinkel,  assistant 
to -Dr.  Helmut  Moritz,  Professor,  Technical  University  at  G  raz  and 
Adjunct  Professor,  Department  of  Geodetic  Science  of  The  Ohio  State 
University,  under  Air  Force  Contract  No.  F  19628-79-C-0075,  The 
Ohio  State  University  Research  Foundation,  Project  No.  711715, 
Project  Supervisor,  Urho  A.  Uotila,  Professor,  Department  of 
Geodetic  Science.  The  contract  covering  this  research  is  admin¬ 
istered  by  the  Air  Force  Geophysics  Laboratory  (AFGL),  Hanscom 
Air  Fo rce  Base,  Massachusetts,  with  Mr.  Bela  Szabo,  Project 
Scientist. 


liooeaeion  Itor 

_J 

mis  GRA*I 

IHDC  TAB 

I1 

.Itoanmounoed 

(Justification _ 

(By 

Distribution/ 

#vai3c_y  *lty  Codc3 


IV 


CONTENTS 

PART  A:  TECHNICAL  BACKGROUND 

Page 

Introduction . 1 

Contouring . 3 

2.  1  Data  sorting  and  retrieving . 3 

2. 2  Prediction  •  . 6 

2.2.1  Inversion-free  prediction. . 6 

2.2.2  Least-squares  prediction . 11 

2.3  Least- squares  regression  . . 16 

2.4  Smooth  surface  representation . 18 

2.4.1  One -dimensional  cubic  spline . 20 

2.4.2  Two-dimensional  cubic  spline . 26 

2.5  The  frequency  content  of  a  bicubic  spline  surface.  .  .  .29 

2.5.  1  Spectrum  of  the  cubic  spline . '  .  .30 

2.  5.  1.  1  Spectrum  of  the  differentiated 

cubic  spline.  . . 40 

2.5.2  Spectrum  of  the  bicubic  spline . .  .  .44 

2.6  Practical  1-D  and  2-D  spline  routines . 46 

2.6.  1  The  cubic  spline  routine  .  . . 46 

2.6.2  The  bicubic  spline  routine . 48 

2.6.3  Interpolation /differentiation  of  splines . 51 

2.7  Contour  finding . 54 

< 

2.8  Optional  contour  procedures . 67 

2.8.1  Contour  smoothing . 67 

2.8.2  Contour  mapping . 70 

2.8.3  Contour  labelling . 71 

2.8.4  Contours  within  a  window . 72 

2.8.5  Contour-free  regions  of  arbitrary  shape . 74 


V 


Page 

2.8.6  Region  boundary  plot* . 85 

2.8.7  Plot  of  horizontal  and  vertical  axes . 88 

2.8.8  Plot  of  a  grid  superimposed  on  the 

contour  plot  . . 89 

2.8.9  Plot  of  data  superimposed  on  the 

contour  plot . .....91 

2. 8!  10  Title  and  label  control  and  plot . .  *92 

2.8.11  Contours  of  surface  derivatives . *95 

3.  Profiles . 97 

3.  1  Explicitly  defined  profiles  . . 98 

3.  1.  1  Optional  profile  procedures . 98 

3. 2  Multiple  profiles.  . . 101 

3. 3  Implicitly  defined  profiles . 101 

3.  3.  1  Implicit  profile  procedures . 104 

4.  Three-dimensional  surface  representations . 105 

4.  1  The  projection . 107 

4. 2  Scale  and  shift . 110 

4,  3  The  3-D  plot . Ill 

4.4  Additional  and  optional  procedures . 112 

PART  B:  ISOLATED  PROBLEMS 

1.  Organization  of  data  . . 115 

2.  Prediction . 119 

3.  Regression.  . . 127 

3.  1  Regression  polynomial  based  on  irregularly 

distributed  data . 128 

3.  2  Regression  polynomial  based  on  regularly 

4 

distributed  data . 130 

4.  Smooth  surface  representations . 133 

4.1  C alculation  of  spline  defining  values . 133 

4.  2  Smooth  surface  interpolation/differentiation . 135 


4.3  Data  approximation  errors 

5.  Axis  plot . 

6.  Grid  plot . 

7.  Title  plot . 

8.  Range  division . 

9.  Clipped  boundary  line  plot  .  • 


Page 
.  141 


.  152 


References  • 


1.  INTRODUCTION 


Technological  developments  during  the  last  decades  have  provided 
us  with  sophisticated  devices  which  are  able  to  collect  data  automatically 
(satellite  altimetry,  Doppler  methods,  satellite  to  satellite  ranging)  or 
semi- automatic  ally  (inertial  navigation)  with  an  enormous  output  rate 
(e.g.,  almost  106  GEOS-3  altimeter  data).  The  tendency  is  clear:  to 
reduce  or  eliminate  manpower  as  much  as  possible  in  the  laborious  and 
expensive  data  acquisition  process.  The  time  seems  to  be  not  far  away 
when  this  first  step  will  have  reached  its  perfection. 

In  the  following  steps,  however,  existing  knowledge  enters  into  the 
data  processing  stage  and  with  this  knowledge  enters  the  human  intellect. 
In  order  to  make  data  processing  faster,  computers  have  been  used 
which  were  designed  in  such  a  way  as  to  permit  the  solution  of  many 
different  kinds  of  problems.  A  relatively  small  number  of  functions  is 
provided  by  the  system.  It  is  usually  up  to  the  user  to  write  his  own  pro¬ 
grams  for  his  own  purposes.  Rigorously  tested  subroutines  are  imple¬ 
mented  into  program  libraries.  A  set  of  programs  might  be  combined 
to  a  large  module  which  would  be  supposedly  capable  of  handling  a 
whole  bunch  of  problems.  Such  kinds  of  modules  are,  in  general  , 
very  large  (e.g.,  a  satellite  orbit  prediction  module)  and  not 
quite  transparent  to  the  user;  therefore,  it  is  highly  desirable 
that  such  a  module  is  intelligent  by  itself.  This  means  that 

the  module  should  be  able  to  check  the  input  (data,  parameters)  for  con¬ 
sistency,  make  adjustments  if  necessary  and  assign  proper  default 
values  to  undefined  parameters.  Such  a  requirement  makes  it 
necessary  to  investigate  the  way  in  which  decisions  are  made  by 
intelligent  human  beings,  to  abstract  this  thinking  process,  and  to 
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tran slate  it  into  a  computer  language.  It  is  most  interesting  and  often 
difficult  to  split  up  fully  automatic  human  decision  processes  into 
steps  or  statements  and  to  find  the  interplay  between  the  visual  per- 
ception  and  the  processing  in  the  brain.  Atypical  example  is  the 
suppressing  of  drawing  contours  within  certain  regions.  We  shall  dis¬ 
cuss  this  problem  in  more  detail  in  Chapter  2.8.5. 

The  very  complexity  of  such  a  module  usually  makes  it  hard  for 
the  user  to  understand  --  this  might  be  the  reason  why  one  speaks 
frequently  about  "black  boxes".  GSPP,  the  Geodetic  Science  Plotting 
Package,  is  such  a  black  box  which  is  designed  for  the  purpose  of 
graphical  representation  of  data  and  smooth  surfaces.  It  is  a  fully 
automatic  link  between  the  stages  of  data  acquisition  and  interpretation 

of  results.  In  view  of  its  great  complexity  and  versatility  it  is  very 
smart  and  simple  to  use;  this  simplicity  is  mainly  due  to  the  control 
part  of  GSPP  which  checks  all  parameters  for  consistency,  makes 
necessary  corrections,  and  assigns  default  values  to  undefined  parameters. 

Since  one  expects  GSPP  to  be  fed  a  large  amount  of  data,  all 
operations  have  been  carefully  checked  and  optimized  in  terms  of  CPU¬ 
time. 

Geodesists  are  used  to  dealing  with  data  on  and  outside  the  surface 
of  the  earth  and  to  preparing  contour  maps  of  surfaces  like  terrain 
contour  maps,  gravity  anomaly  contour  maps  or  geoid  contour  maps. 

These  two-dimensional  representations  of  surfaces  have  --  besides 
providing  a  general  behavior  of  the  surface  --  the  advantage  of  allowing 
the  user  to  also  interpolate,  to  some  extent,  the  information  contained  in 

the  contours.  This  makes  a  contour  map  superior  to  a  single  pro¬ 
file  as  far  as  global  numerical  information  is  concerned  and  also 
superior  to  a  three-dimensional  view  which  provides  the  user  a  unique 
spatial  image  but  lacks  numerical  information. 

Therefore,  contour  maps  are  indeed  unique  and  this  is  the  reason 
why  the  present  report  starts  with  the  discussion  of  contour  maps  and 
not  with  profiles  as  one  would  expect. 
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The  user  primarily  interested  in  applications  may  skip  the  technical 
Part  A  which  is  intended  to  throw  light  into  the  darkness  of  the  black 
box;  he  may  start  with  Part  B  and  consult  Part  A  when  the  need  arises. 

2.  CONTOURING 

Before  we  start  with  a  detailed  description  of  the  whole  contouring 
process  it  should  be  pointed  out  that  GSPP  is  not  designed  for  the  purpose 
of  representing  terrain  structures  with  all  its  many  details  as  a  digital 
terrain  model  will  do  --  it  is  designed  for  representing  smooth  surfaces 
without  artificial  structures.  Such  smooth  surfaces  are  derived 
from  surface  data  located  at  the  grid  points  of  a  regular  rectangular 
grid.  In  almost  all  practical  applications  we  are,  however,  far  away 
from  this  ideal  situation,  for  three  reasons:  first,  the  surface  information 

is  sampled  at  points  which  are,  in  general,  irregularly  distributed 
(terrain  heights,  gravity  anomalies,  geoidal  heights);  second,  the  data 
are  often  heterogeneous  in  nature  and  contain  surface  information  only 
implicitly  (data  combinations  in  physical  geodesy);  third,  the  data  are 
usually  disturbed  by  some  kind  of  noise.  All  these  deficiencies  make 
contouring  a  non-trivial  and  also  non-unique  task. 


2.  1  Data  sorting  and  retrieving 

Theoretically,  all  available  data  should  be  used  in 
order  to  really  achieve  a  prediction  with  the  minimum 
variance.  Practically,  this  is  neither  possible  because  of  the  enormous 
amount  of  geodetic  data  collected  so  far,  nor  is  it  necessary:  a  pre¬ 
dicted  value  at  a  point  P  depends  primarily  on  the  data  in  the  neighbor¬ 
hood;  remote  data  will  often  contribute  very  little  to  the  result. 

Therefore,  one  accepts  the  sub-optimal  solution  and  takes  only  a  relatively 
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small  number  of  data  for  a  single  prediction  into  account.  All  other 
data  which  exceed  a  certain  distance  from  the  prediction  point,  are  not 
considered. 

This  would  make  the  calculation  of  all  distances  between  all  data 
points  and  all  prediction  points  necessary,  a  very  time  consuming  task 
especially  when  working  on  the  surface  of  the  sphere  where  "expensive" 
trigonometric  functions  are  necessary  in  order  to  calculate  the  distance. 
For  this  reason,  it  is  absolutely  necessary  to  order  the  data  according 
to  some  pre-  selected  pattern  in  order  to  be  able  to  single  out  the 
unnecessary  part  in  a  simple  and  fast  way.  This  is  also  one  of 

the  properties  of  data  bases.  It  is  not  our  intention  to  establish  a 
sophisticated  data  base  with  a  very  complicated  tree  structure;  all  we 
want  to  do  is  to  find  a  simple  way  of  ordering  some  irregularly  distributed 
data. 

For  the  sake  of  simplicity,  we  assume  data  to  be  1-point  data. 

Let  (xj,  yj),  i  =  1 . n  be  the  coordinates  of  a  data  point  P.  in 

a  Cartesian  coordinate  system.  Then  an  obvious  way  of  arranging  the 
data  would  be  to  select  a  certain  "working  domain",  which  is  prefer¬ 
ably  a  rectangle,  by  defining  the  boundaries  of  the  rectangle  (x,  ,  x  , 

U 

y_  ,  y  ),  the  lower  and  upper  x-  and  y-coordinates  if  the  rectangle 
Li  U 

happens  to  be  parallel  to  the  coordinate  lines.  This  rectangle  will 
then  be  subdivided  into  a  number  of  subrectangles.  The  idea  of  the  data 
sorting  process  is  to 

1.  findforeach  data  point  the  corresponding  subrectangle  (element), 

2.  count  the  number  of  data  within  each  element,  and 

3.  generate  pointer  vectors  which  allow  us  to  find  all 
data  points  within  a  specific  element. 

Note,  that  all  data  remain  at  their  original  storage  locations  --  only 
additional  information  is  produced. 

Let  us  now  follow  these  three  steps  in  detail: 
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Assume  the  working  domain  is  subdivided  into  I  #  J  elements 
(i  inx-,  J  in  y-direction)  and  enumerate  these  elements  such  that 
an  integer  (i-l)*J  +j  is  assigned  to  the  element  (i,  j).  Then  for 
each  data  point  there  exists  an  integer  which  identifies  an  element  if  all 
data  are  located  within  the  working  domain;  this  vector  will  be  denoted 
by  a;  the  value  of  a^.  depends  only  on  the  position  of  P^, 

\  ■  £<pk)  • 


The  length  of  a  is  equal  to  the  number  of  data  points.  At  the  same 
time  a  counter  vector  b  (length  =  number  of  elements)  counts  the 
number  of  data  points  found  in  each  element.  After  all  data  have  been 
located,  b  contains  information  about  the  total  number  of  data  within 
each  element. 

The  problem  is  to  generate  pointer  vectors  such  that  it  is  possible 
to  find  all  data  within  a  prescribed  element.  For  this  purpose  another 
vector  c  has  to  be  generated  which  has  the  same  length  as  b  and  is  its 
partial  sum.  It  provides  the  information  of  the  starting  place  in  the  final 
pointer  vector  d  such  that  the  first  data  point  in  a  particular  element 
Z  =  (i-  1)  #  J  +  j  has  the  index 

m  =  d[c(4)]  , 

the  second  data  point  in  the  same  element  l  has  the  index 
m  =  d  [c(-2)  +  1  ]  , 

the  last  data  point  in  element  Z  has  the  index 
m  =  d  [c4)  +  b(X)  -  l]  . 

If  the  working  domain  does  not  enclose  all  the  data,  the  vectors 
b  and  c  have  to  be  longer  by  1  element  such  that  the  number  of  data 
outside  the  working  domain  can  be  stored.  In  this  case  b(l#J  +1) 
stores  the  number  of  data  outside  the  working  domain  and  d[c(l*J  +1)] 
is  the  index  of  the  first  data  point  found  outside  the  working  domain. 
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The  program  responsible  for  this  data  organization  job  is  called 
OAF;  it  is  independent  and  can  be  run  separately.  The  time  for  the 
organization  process  depends  almost  entirely  cn the  number  of  data  and 
increases  linearly  with  the  number  of  data.  In  order  to  give  an  idea 
about  the  efficiency  we  give  one  example:  to  order  100  000  data 
according  to  the  process  described  above  takes  1.  5  sec  CPU-time  on 
a  AMDAHL  470V/6-II  computer  which,  in  the  opinion  of  the  author, 
is  quite  fast. 

2.  2  Prediction 

In  general,  most  kinds  of  data  will  be  irregularly  distributed.  For 
many  reasons,  however,  a  regular  distribution  (preferably  on  a  regular 
rectangular  grid)  would  be  welcomed: 

a)  regularly  distributed  data  admit  a  simple  data  manage¬ 
ment  and  an  easy  data  retrieving; 

b)  a  regular  data  distribution  (on  a  rectangular  grid)  is 
particularly  important  if  one  wants  to  fit  a  smooth  surface 
like  a  bicubic  spline  function  to  the  data; 

c)  a  regular  distribution  is  necessary  for  fast  post-process*' 
ing  of  data  using  techniques  like  the  Fast  Fourier  Transform 
method;  but  also  the  standard  least- squares  methods 
would  gain  considerably  in  terms  of  computer  time  by 
taking  symmetries  of  the  covariance  matrix  into  account. 

There  are  a  number  of  possibilities  of  obtaining  predicted  (or  inter¬ 
polated)  values.  We  will  be  discussing  only  two  of  them,  the  inversion- 
free  prediction  and  least- squares  collocation. 

2.2.1  Inversion-free  prediction.  In  many  cases  we  face  the 
following  situation:  there  is  given  a  huge  number  of  data,  irregularly 
distributed,  of  a  homogeneous  type  and  almost  free  of  errors;  predict  a 
reasonable  smooth  surface.  If  one  would  like  to  obtain  an  optimal 
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solution  in  the  sense  of  least  norm,  a  least- squares  solution  is  appro¬ 
priate.  Minimizing  the  mean  square  prediction  error  leads  to  least- 
squares  prediction  and  least- squares  collocation.  This  method,  however, 
requires  an  a  priori  covariance  function  to  be  known  which  is  missing 
in  many  kinds  of  problems.  Moreover,  an  inversion  of  a  large  matrix 
or  many  inversions  of  small  matrices  make  solutions  of  such  kind 
of  problems  too  often  a  very  expensive  task. 

A  simple  and  cheap  (however,  not  the  best)  alternative  is  to  define 
the  predicted  value  as  a  weighted  average  of  all  data  in  the  neighborhood 
of  the  prediction  point.  The  weights  will  depend  on  the  distances  between 
data  points  and  the  prediction  point.  One  kind  of  such  a  prediction  is 
the  following: 

Let  {f.},  i=  l,...,n  be  n  homogeneous  data;  then  the  predicted 
value  at  a  point  P  is  given  by 


with  s^p  ...  distance  between  data  point  and 
prediction  point  P, 

q  ...  power  of  prediction. 

This  kind  of  prediction  has  the  advantage  of  exactly  reproducing  the 
data  and  of  not  involving  the  calculation  of  time  consuming  functions 
and  matrix  inversions.  This  makes  it  particularly  useful  for  large- 
scale  applications  of  the  type  of  problems  discussed  above. 

It  is  quite  instructive  to  compare  this  kind  of  prediction  with 
least- squares  prediction: 


Assume  we  are  given  3  error  free  geoidal  heights  and  want  to  pre¬ 
dict  a  geoidal  height  at  a  point  P.  Then  the  inversion-free  prediction 
gives  the  result 

£P  =  7“  f*  +  ai  £z  +  *3*3) 

with  0!i ,  i  =  1,2,3  are  weights  depending  only  on  the  distance  between 
P  and  all  P^  (and  not  on  the  distance  between  P^  and  P^);  c  is 

the  sum  of  all  0^  .  In  the  case  f  =  ft  =  f,  =  f3  we  obtain 

fp  =  '  “  *  («  i  +  «2  +  a  3)  =  f  “  =  £ 

and  the  predicted  value  £  is  independent  of  the  position  of  the  prediction 
point;  consequently,  the  surface  generated  will  be  a  horizontal  plane.  In 
the  case  of  a  general  vector  {  f^  }  ,  the  surface  will  not  be  a  plane 
anymore,  however,  its  trend  outside  the  data  region  (here  a  triangle) 
will  be  to  be  a  plane  with  function  value  equal  to  the  mean  value  of  all 
data.  The  reason  is  that  the  differences  between  the  weights 

6<*  -  j;  <tt)  «s  -  -9  — « 

6  s 

tend  to  zero  because  of  the  factor  —  . 

s 

The  least- squares  prediction  behaves  essentially  different 
(Cij  ...  covariances), 
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Under  the  assumptions  made  above  (  ffj]  =  error-free  geoidal 
heights),  the  main  diagonal  elements  will  be  constant  and  equal  to  the 
variance  C0.  Assume  now  for  a  moment  that  the  correlation  length  of 
the  covariance  function  C  becomes  very  small  compared  with  the 
smallest  initial  distance  between  the  data;  then  the  covariance  matrix 


will  be  highly  dominated  by  the  diagonal  and  its  inverse  will  only  slightly 


deviate  from  a  purely  diagonal  matrix  with  diagonal  elements  being 

almost  equal  to  the  inverse  of  the  variance,  -7-  .  The  covariances 

—  C° 

Cp^,  Cp^i  Cp3  °e  small  *°r  large  FF.»  as  a  consequence  the 


surface  will  become  wavy  in  between  the  data  points  and  tend  to  zero 
outside  the  data  region,  (In  the  spherical  case  the  surface  will  be 
such  that  its  integral  over  the  sphere  vanishes.  )  For  larger  correlation 
lengths  the  predicted  surface  will  become  stiff er.  In  any  case,  the 

covariance  function  controls  the  surface;  the  "reproduction"  of  the 
surface  on  the  basis  of  its  sampled  data  consists  basically  of  its 
"inversion". 

Therefore,  one  main  difference  between  both  methods  is  that  in 


least- squares  prediction  the  surface  is  generated  by  actually  inverting 
the  surface  information;  in  the  inversion-free  surface  prediction, 
correlations  between  the  data  are  completely  neglected. 

A  very  interesting  light  has  been  thrown  recently  on  the  inversion- 
free  prediction  (interpolation)  in(5unkel,  1980).  It  has  been  shown  that  (2-1) 


can  be  represented  in  terms  of  base  functions  B  (x), 


n 

L  f,  Bq  (P-Pj)  . 


Each  base  function  B  (x)  can  be  shown  to  be  identical  with  the  Fourier 
transform  of  a  polynomial  spline  of  degree  q-1.  Since  the  Fourier 
transform  of  a  polynomial  spline  approaches  a  rectangle  with  increasing 
degree  q,  the  base  functions  also  approach  a  rectangle  as  the  power 
of  prediction  increases.  If  q  ,  the  "interpolation"  function  defined 
by  (2.  1)  will  be  a  pure  step  function. 

Figure  2.  1  gives  an  impression  of  what  happens  if  the  power  of 
prediction  is  much  too  large.  The  data  indicated  by  "+"  are 
"interpolated"  by  functions  which  are  very  close  to  step  functions 
(Adiguzel,  1979). 


FIG.  2.  1  --  Step  function  like  interpolation  for  very  high 

powers  of  prediction. 

On  the  other  hand,  if  the  power  of  prediction  q  is  too  small,  0  <  q  £  1, 
the  interpolation  function  tries  to  average  out  all  data  and,  at  the  same 
time,  it  reproduces  the  data;  this  causes  the  function  to  produce  cusps 
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in  the  neighborhood  of  the  data  points  which  should  obviously  be  avoided. 
Since  neither  steps  nor  peaks  are  what  one  considers  as  a 
reasonable  interpolation  function,  one  of  the  often  accepted  compromises 
is  the  choice  q  =  2  {Schumaker,  1976,  Bybee  and  Beaross,  1978);  other 
users  rather  prefer  a  value  around  3.5  (Bjerhammar,  1973;  Davenport 
--private  communication).  GSPP  allows  the  user  to  make  his  own 
choice  between  the  limit  values  q  =1.5  and  q  =  4.  The  method  de¬ 
scribed  above  goes  back  to  (Shepard,  1964)  and  (Bjerhammar,  1973). 

Another  compromise  is  a  splitting  of  the  weight  function  into  three 
or  more  parts,  each  of  which  shows  a  different  behavior  (has  different 
powers  q  assigned),  „ 


r-i  , 


0<  r  £  | 


(  0  R  <  r 

This  weight  function,  proposed  by  (Shepard,  1964),  is  continuous  and 
continuously  differentiable  and  vanishes  outside  r>  R,  the  radius  of 
prediction.  Consequently,  the  choice  of  R  controls  the  interpolation 
behavior  and,  at  the  same  time,  defines  the  prediction  circle.  (All 
data  within  this  circle  of  radius  R  around  the  prediction  point  contribute 
to  the  predicted  value.  )  This  function  is  the  default  function  in  GSPP 
if  no  power  of  prediction  is  defined  by  the  user. 

2.2.2  Least- squares  prediction.  The  inversion-free  prediction 
as  discussed  in  the  foregoing  chapter  is  applicable  in  a  very  special 
case,  when  data  are  homogeneous  and  free  of  noise.  These  assumptions 
are  often  almost  satisfied,  however,  the  general  case  of  heterogeneous 
noisy  data  cannot  be  treated  with  inversion-free  prediction. 

As  stated  in  the  introduction,  the  module  described  here  is  primarily 
intended  for  applications  in  geodesy,  more  specifically,  for  purposes 
of  the  determination  of  the  gravity  field  of  the  earth.  There  are  many 
types  of  data  used  in  physical  geodesy;  all  of  them  bearing  information 


I 


about  the  gravity  field,  all  of  them  are  more  or  less  noisy.  The  problem 
is  to  optimally  combine  these  data  in  such  a  way  that  the  gravity  field 
determined  on  the  basis  of  these  data  deviates  from  the  true  one  as 
little  as  possible.  Least- squares  collocation  turned  out  to  be  particularly 
useful  for  the  solution  of  such  kind  of  problems.  The  goodness  of 
the  predicted  field  depends  considerably  on  the  covariance  function  intro¬ 
duced  which  should  match  the  average  behavior  of  the  gravity  field  as 
close  as  possible  (Schwarz  and  Lachapelle,  1979).  For  reasons  of 
continuity,  the  covariance  function  will  be  briefly  described. 

The  general  form  of  a  homogeneous  and  isotropic  covariance 
function  of  the  disturbing  potential  can  be  expressed  by 

Y  /rb\  n+1 

K(P,Q)  =  L  ^  \rr' )  Pjcos  *  )  (2-2) 

n=N0 

with  P,Q  =  points  outside  the  sphere  r  =  Rg  , 

r,  r'  =  geocentric  radii  of  P  and  Q, 

0  pQ  =  spherical  distance  between  P  and  Q, 

=  positive  coefficients, 

Pn(cos  0)  =  Legendre  polynomial  of  degree  n, 

N0  =  starting  value  of  the  summation  (N0  ^  2), 

Rg  =  radius  of  the  Bjerhammar  sphere. 

K(P,Q)  is  symmetric  with  respect  to  P  and  Q  and  harmonic  outside 
the  sphere  r  =  Rg. 

Geodetic  measurements  are,  in  general,  nonlinear  functionals  of 
the  gravity  field  (-t  station  position).  After  linearization  they  can  be 
expressed  as  linear  functionals  of  the  anomalous  gravity  field  (+  station 
coordinates). 

In  order  to  predict  a  linear  functional  of  the  gravity  field  at  an 
arbitrary  point,  it  is  first  of  all  necessary  to  establish  the  covariance 
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matrix  which  represents  so  to  say  the  structural  relation  between 
the  gravity  field  and  the  data.  An  element  of  the  covariance  matrix 
C„  is  the  result  of  an  operation  which  maps  the  covariance  function 
into  the  real  number  line;  the  mapping  consists  of  applying  the  linear 
operations  and  corresponding  to  the  data  i  and  j  on  the  covariance 
function  K  : 


C..  =  L.L.K. 

ij  i  J 

In  the  same  way  are  the  covariances  between  the  predicted  quantity 
and  the  data  obtained: 


'Pi 


*  W- 


Denoting  the  (linearized)  vector  of  measurements  by  f  as  above,  one 

can  find  a  simple  linear  relation  between  predicted  quantities  and  data 

T  - 1 

(Moritz,  1978):  fp  =  Cp  C"  f  ,  or  more  detailed, 


fP  =  tCPl’  Cp2*  ’  *  *  ’  CpJ 

Cu 
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Taking  also  noise  into  account  one  has  to  add  to  C  the  corresponding 
error  covariance  matrix  of  the  data  which  usually  happens  to  be 
diagonal.  In  the  general  case  with  incorporated  model  parameter 
determination,  formula  (2-3)  changes  slightly;  however,  we  shall 
limit  ourselves  to  the  case  discussed  above.  The  error  covariances 
between  the  predicted  quantities  at  the  points  P  and  Q  can  be  shown 
to  equal 

“PQ  =  CPQ  -  CPC‘‘CQ  ! 

the  variance  is  obtained  for  Q  =  P  (Heiskanen  and  Moritz,  1967;  P  269  ff-. ) 
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The  covariances  can  be  derived  from  a  model  covariance 
function  which  should  be  simple  in  order  to  keep  computation  time  low; 
on  the  other  hand  it  should  have  certain  properties  (represented  by- 
essential  parameters)  which  are  determined  by  the  general  features  of  the 
gravity  field.  These  two  requirements  exclude  each  other.  Therefore,  it  is 
generally  a  rather  time  consuming  task  to  find  numerical  values  for 
all  covariances  involved,  even  when  using  well-designed  algorithms 
(Tscherning,  1976). 

Therefore,  it  has  been  investigated  if  one  could  get  rid  of  this  burden 
by  taking  advantage  of  finite  elements  and  approximate  the  covariance 
function  (Sunk el,  1978).  The  idea  behind  this  is  simple  and  has  been 
frequently  applied  in  many  fields,  the  network  principle:  generate  a 
net  of  fixed  points  (here  grid  points  in  a  two-dimensional  space)  and 
perform  very  accurate  measurements  at  these  points  (here,  calculate 
exact  covariances);  these  fixed  points  serve  as  a  basis  for  small  scale 
measurements  which  can  be  obtained  from  simple  devices  (here,  differentia¬ 
tion  -  interpolation  of  finite  elements  representing  the  covariance  func¬ 
tion  in  a  certain  range).  The  bicubic  spline  function  turned  out  to  be 
particularly  useful  for  such  a  purpose.  In  (Siinkel,  1979  )  such  a 
covariance  approximation  procedure  has  been  described;  a  FORTRAN  IV 
computer  program  which  has  been  designed  for  such  purposes,  is  used 
by  GSPP.  This  program  allows  us  to  obtain  covariances  of  second  and 
lower  order  derivatives  of  the  disturbing  potential  at  a  CPU-time  level 
of  some  3*  10 ~4  sec  on  a  IBM  370  system.  This  advantage,  however, 
is  to  some  extent  balanced  by  disadvantages  of  one  needing:  to  know  the 
maximum  range  (in  space)  in  which  one  is  working,  to  generate  a 
network  of  covariances  and  store  this  network  on  a  (advantageously 
permanent)  file  before  one  can  actually  call  the  covariance  approximation  sub¬ 
routine;  and  finally  to  realize  covariances  differ  slightly  from  the  exact  ones 
derived  from  the  model  covariance  function;  the  differences,  however, 
can  be  kept  arbitrarily  small.  This  is  the  price  we  pay  for  a  90  ^ 


reduction  in  the  computation  time  for  covariances. 

It  should  be  clear  that  such  sophisticated  procedures  require 
decision  making  before  actually  using  the  module  GSPP.  This  means 
that  the  covariance  approximation  procedure  should  be  used  only  for 
large  scale  applications  --  for  such  a  purpose  it  was  designed.  In 
such  cases  it  really  pays  back  by  saving  enormous  amounts  of  computer 
time.  For  medium  and  small  scale  applications  one  can  use  the  exact 
covariances  derived  from  some  model. 

Detailed  information  about  the  covariance  approximation  algorithm 
can  be  found  in  (Stinkel,  1979  )s  the  principle  can  be  briefly  described 
in  the  following: 

An  important  property  of  the  covariance  function  which  makes  a 

two-dimensional  approximation  for  theoretically  all  points  outside  the 

Bjerhammar  sphere  possible,  is  its  dependence  on  essentially  two  variables, 

the  spherical  distance  $  and  the  product  rr'.  Since  cos0  can  vary  only 

between  -1  and  +1  and  Ri  /  rr1  has  a  minimum  value  of  0  for  r  * 

B 

and  a  maximum  of  1  for  r  =  r1  =  Rg,  the  covariance  function's  domain  of 
definition  is  the  rectangle 

[  - 1  =  fc  —  1  ,  0<s<l] 


with 


t  =  cos  0 


ik 

rr' 


Since  practically  all  geodetic  operations  are  performed  on  or 
close  to  the  surface  of  the  earth,  the  domain  of  definition  is  reduced 
considerably  for  all  practical  applications,  e.  g.,  working  within  a 
spherical  distance  range  of  0  s  ^  s  10°  and  within  an  altitude  range 
from  0  to  300  km,  the  domain  reduces  to 


[  0.985  *  t  s  1  »  0. 999  >  s  >0.912]  . 

Once  the  user  has  made  his  decision  about  the  r'anges  in  which  he 
is  going  to  work,  a  rectangular  grid  in  s  and  t  can  be  arranged.  A 
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specific  program  has  been  designed  for  the  calculation  of  different  kinds 
of  covariances  at  all  grid  points.  This  network  of  covariances  has 
to  be  calculated.  The  calculation  itself  is  a  three  stage  process:  in 
the  first  part  the  grid  element  corresponding  to  the  particular  data 
location  has  to  be  found;  in  the  second  stage  the  necessary  interpolations 
and  differentiations  of  the  covariance  network  are  performed;  and  in 
the  last  part  the  program  calculates  the  partial  derivatives  of  cos  $ 
with  respect  to  the  spherical  coordinates  at  the  particular 

points  P  and  Q.  So  far  as  covariance  approximation  is  concerned. 

The  prediction  program  itself  has  control  over  the  prediction 
region,  a  rectangle  surrounding  the  prediction  point.  This  region  can 
be  defined  by  the  user  in  terms  of  a  radius  of  prediction.  As  long  as 
the  program  finds  a  minimum  number  of  data  for  a  single  prediction,  it 
is  fine;  if  this  is  not  the  case  the  prediction  region  will  be  enlarged 
until  a  minimum  number  of  data  have  been  found.  Three  different 
kinds  of  prediction  are  possible,  the  inversion-free  prediction  for  homo¬ 
geneous  and  error-free  data,  least- squares  prediction  with  accurate 
covariances  derived  from  a  model  covariance  function,  and  least- 
squares  prediction  with  approximated  covariances.  In  the  case  of  least- 
squares  prediction,  the  system  returns,  apart  from  the  predicted  grid, 
also  the  root  mean  square  prediction  errors.  The  calculation  of  these 
errors  can  be  suppressed.  The  predicted  grid  represents  the  surface 
insofar  as  the  bicubic  spline  interpolation  function  based  on  this  grid 
is  interpreted  as  the  surface  in  consideration. 

2.  3  Least- squares  regression 

A  surface  predicted  by  one  of  the  methods  discussed  above  is 
capable  of  Representing  even  small  and  local  details  of  the  surface  if  the 
data  contain  such  an  information.  Such  detailed  representations  can 
only  be  described  by  a  large  number  of  parameters;  in  case  of  least 
squares  prediction  the  number  of  parameters  (coefficient^  is  equal  to 
the  number  of  data  .  For  many  reasons,  however,  one  is  often  not 
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only  interested  in  the  local  details  but  also  in  the  global  features  of 
the  field  represented  by  the  data.  Such  global  features  can  be  described 
by  a  relatively  small  number  of  parameters.  The  problem  consists  of 
the  selection  of  the  model,  which  is  usually  a  polynomial  of  some  degree, 
in  establishing  the  relation  between  the  data  and  the  model  parameters 
and  in  the  solution  of  the  linear  system  of  a  size  equal  to  the  number  of 
parameters.  The  most  obvious  solution  is  a  least- squares  solution, 
in  this  context  we  speak  about  least- squares  polynomial  approximations 
or  simply  about  least- squares  regression.  Typical  examples  of 
practical  applications  are  trend  eliminations  in  all  natural  sciences 
(e.g.,  determination  of  spherical  harnomic  coefficients  up  to  and  including 
degree  and  order  N  based  on  data  sensitive  to  the  earth's  gravity  field). 

The  basic  idea  comes  from  an  old  interpolation  theorem,  the 
century-old  Weierstrass  approximation  theorem  and  the  least- squares 
principle.  Very  loo  sly  speaking  the  first  theorem  says  that  all  n 
coefficients  of  a  polynomial  model  can  uniquely  be  determined  from  n 
independent  data  linearly  related  to  the  model.  Weierstrass'  theorem 
asserts  that  a  continuous  function  can  be  uniformly  approximated  by  a 
polynomial  on  a  closed  interval.  The  least- squares  principle  guarantees 
the  uniqueness  and  existence  of  a  shortest  distance  between  a  point 
(vector  of  data)  and  the  hyperplane  spanned  by  the  linearly  independent 
base  functions  (polynomials  1,  x,  y,  x2,  xy,  y2,  . . . ,  y  ).  The 
parameters  {a.}  can  be  immediately  found  from  the  solution  of  the 
normal  equations: 

<  F  -  $  a  ,  $  >  =  0  (2. 3- 1) 

with  F  =  { f  ^  }  =  [Li.f  }  data  (L.  ..linear  functional) 

a  =  f  a^  }  . . .  polynomial  coefficients 

<t>  =  {  }  =  {  Li  v?  j }  ...  design  matrix 

(  ¥>j  ...  base  functions  ). 


*  (Davis,  1975,  pp.  24  ff. ,  107  ff . ,  158  ff . ) 
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"With  an  a  priori  error  covariance  matrix  P^of  the  data,  the  well 
known  solution  for  the  parameters  is 

a  *  (  $TP  $)  1  $TPF  ;  (2.3-?.) 

the  individual  data  reproduction  errors  are  given  by 

AF  =  F  -  fa.  (2.3-3) 

E  the  number  of  parameters  is  equal  to  the  number  of  data,  AF 
vanishes  identically  (simple  interpolation). 

El  order  to  avoid  any  misinterpretation,  it  should  be  pointed  out 
that  least-squares  polynomial  approximations  do  by  no  means  replace  or 
compete  with  prediction  solutions  based  on  the  least- squares  collocation 
principle  or  any  other  prediction  method;  they  supplement  these  solutions 
insofar  as  they  provide  a  trend  information.  Such  trend  calculations  may 
be  quite  useful  for  a  number  of  problems;  however, 
the  user  should  be  warned  not  to  work  with  a  high  degree  trend  poly¬ 
nomial  and  to  make  sophisticated  interpretations  on  the  basis  of  the 
results:  polynomials  show  the  tendency  to  oscillate  between  data  and 
do  not  hesitate  to  show  completely  abnormal  features  in  data  free  regions, 
(see  as  an  example  Fig.  2.2). 


2.4  Smooth  surface  representation 

To  predict  values  at  points  where  no  measurement  has  been  per¬ 
formed,  is  usually  very  expensive  in  terms  of  computer  time.  The 
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reasons  are  multifold:  it  is  necessary  to  compute  distances  between 
each  calculation  point  and  theoretically  all  data  points  (this  deficiency 
is  shared  by  all  prediction  methods);  the  use  of  least- squares  collocation 
usually  involves  the  calculation  of  trigonometric  and/or  logarithmic 
functions,  and  the  multiplication  of  large  vectors  with  large  matrices 
for  each  single  prediction.  (We  don't  consider  at  this  point  the  inversion 
of  the  covariance  matrix  because  this  is  independent  on  the  number  of 
predictions.)  Therefore,  it  is  absolutely  necessary  to  find  a  representa¬ 
tion  of  the  surface  which  is 

a)  based  on  the  data  and/or  the  predicted  values, 

b)  smooth, 

c)  local, 

d)  simple. 

Requirement  (a)  is  evident;  a  sufficient  degree  of  smoothness  is  desired 
in  order  to  admit  surface  differentiations  (slope  maps,  etc.), 
the  interpolating  function  should  be  sensitive  to  a  point  disturbance  only 
in  its  neighborhood  which  is  referred  to  as  a  local  behavior;  last, 
but  by  no  means  least,  should  the  interpolating  function  be  simple  in 
order  to  make  the  interpolation  /differentiation  process  fast. 

Naturally,  there  does  not  exist  a  function  which  fulfills  all  these 
requirements  fully.  Single  polynomials  are  not  local,  linear  interpolating 
elements  are  not  smooth.  An  optimal  compromise  is  possibly  a  bicubic 
spline  function  which  is  sufficiently  smooth  (continuous  second  order 
derivatives),  is  strictly  local,  and  still  a  relatively  simple  interpolating 
element.  A  disadvantage,  however,  is  that  bicubic  spline  functions  are  based 
on  a  regular  rectangular  grid;  smooth  surfaces  based  on  irregularly  distrib 
uted  data  are  possible;  the  computational  effort,  however,  is  huge.  On  the 
other  hand,  a  regular  rectangular  grid  of  data  is  anyway  the  most  natura’  way 
of  storing  two-dimensional  information,  and  therefore,  this  restriction  loo ses 
much  of  its  power.  In  the  following  we  give  for  the  sake  of  completeness  a 
short  description  of  spline  interpolation  in  one  and  two  dimensions  which  is  a 
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short  outline  of  (Sunk el,  1980)  and  the  procedure  GSPP  uses  to  smoothly 
represent  a  surface  based  on  gridded  data. 

2.4.  1  One-dimensional  cubic  spline. 

The  one -dimensional  cubic  spline  is  a  basic  tool  not  only  for 
representing  a  smooth  curve,  but  also  for  surface  fitting.  The  cubic 
spline  is  just  one  of  an  infinite  number  of  splines;  its  pleasant  proper¬ 
ties  make  it  unique  among  all  splines.  It  is  a  function  which  is  twice 
continuously  differentiable  and  therefore  very  smooth.  Basic  for  the 
cubic  spline  is  the  cubic  basis- spline,  or  simply  B- spline,  defined  on  a 
grid  with  constant  grid  spacing  equal  to  1.  (in  such  a  case  one  usually 
speaks  of  cardinal  splines  because  of  its  definition  on  the  sequence  of 
cardinal  numbers.)  Such  a  cardinal  B-spline  is  a  piecewise  cubic 
polynomial  with  bounded  support;  it  is  twice  continuously  differentiable 
on  the  whole  real  line  (-«.  •»  ).  Centered  at  the  zero  point,  it  can  be 
expressed  by 


t 


B(x  ) 


7  2.  (  -l)k  (  J)  (x  +  2  -  k)J 

6  k=0 


(2.4-1) 


with 


x  for  x  >  0 
0  otherwise 


(I.  J.  Schoenberg,  1973,  p.  11).  Its  support  is  the  open  interval 
(-2,2).  Explicitly  written,  B(x)  satisfies  the  equations 

13 1 x J 3  -  6x2  +4  for  Osjxj^  1 

-|x|3  +  6x2  -  12|xj+8  for  l£|x|^  2 

0  2£|x| 

It  can  be  seen  immediately  that  B(x )  is  symmetric:  B(x)  =  B(-x)  . 
Its  function  values  at  the  knots  are 


B(±2)  =  0,  B(±l)  =  |,  B(0)  =  | 


(2.4-2) 


Figure  2.3  shows  the  cubic  cardinal  B-spline  together  with  its 
derivatives  up  to  and  including  order  3  which  is  a  step  function. 
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Given  function  values  f,  and  second  derivatives  f. "  at  all 

k  k 

integers  on  the  real  line,  it  can  easily  be  shown  that 


F  (x)  =  A  cfc  B  (x-k) 


(2.4-3) 


reproduces  the  data  and  is  a  unique  interpolating  function  (SUnkel,  1980). 
The  coefficients  {c  ^5  can,  in  fact,  be  expressed  by  c  ^  =  f^-  Af  it  , 

In  almost  all  practical  applications,  however,  the  second  derivatives 
are  not  known.  Therefore,  the  question  arises  whether  it  is  possible 
to  find  an  interpolating  cubic  spline  which  is  only  defined  on  the 
function  values  (f  }  at  the  grid  points.  Or,  formulated  differently, 

fcC 

one  would  like  to  have  functions  S(x)  such,  that 


F(x)  =  .4  fkS(x"k) 


(2.4-4) 


is  a  cubic  interpolating  spline.  Such  a  function  does  exist 

(i.  J.  Schoenberg,  1973);  it  can  be  expressed  as  a  discrete  convolution 


of  the  form 


(x)  =  Z  a. 


(2.  4-5a) 


the  coefficients  {<7. }  can  be  determined  from  the  condition  that 

J 


sw  ■  \  -{in;?. 


(2. 4-5b) 


Since  the  cubic  B-spline  has  non-vanishing  function  values  at  only  3 
knots  (-1,0,  1),  it  follows  with  (2.4-5a,  b)  that  the  infinite  sum  reduces 
to  a  sum  over  only  3  B- splines  for  S(x)  =  S(k): 

S(k)=a  B(  1)  +  cr  B(0)  +  Vk* 


or  explicitly  , 
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S(-2)  = 

or_3B(l) 

+ 

<7-2B(0) 

+ 

tf-i 

b(-  1) 

=  0 

s(-i)  = 

ff.2B(l) 

+ 

*-»B(0) 

+ 

<7o 

B(-  1) 

=  0 

S(  0)  = 

<7-1  B(  1) 

+ 

B(0) 

+ 

<7i 

b(-  1) 

=  1 

S(l)  = 

Q 

o 

td 

»— * 

+ 

al  B(0) 

+ 

^  2 

b(-  1) 

=  0 

S(2)  = 

CTl  B(  1) 

+ 

ff2  B(0) 

+ 

^3 

b(-  1) 

=  0 

Therefore,  in  order  to  find  the  infinite  vector  of  coefficients 

{ (T.  }  .  j  =  . .  • ,  **,  we  have  to  solve  the  infinite  system  of 

3 

equations  (2.4-6)  which,  with  B(k)  from  (2.4-2),  is  given  by 


***** 

•  •  •  •  • 

*  *  *  *  • 

...0  1  4  1  0... 

...  0  1  4  1  0 . . . 

...0  1  4  1  0... 

. .  .0  1  4  1  0. . . 

...  0  14  10... 


• 

• 

• 

• 

• 

CT-2 

• 

0 

0 

<7o 

= 

6 

<7l 

0 

!2 

0 

• 

• 

• 

L  -J 

(2.4-6)' 


The  transformation  matrix  is  of  infinite  dimension,  is  symmetric,  and 
because  of  its  "row- shift"  structure  of  Toeplitz  form,  it  is 
circulant.  Such  kind  of  matrices  are  well  known  to  have  inverses  of 
the  same  type  (R.M.  Gray,  1971);  with  the  help  of  Fourier  techniques 
it  is  fairly  easy  to  find  the  solution  vector  qp  .  In  (H.  Sttnkel,  1980) 
it  is  shown  in  detail  how  the  actual  solution  can  be  obtained;  the 
result  is 


a.  =  /3  (-2  +  /3)  ^ 


(2.4-7) 


Consequently,  the  cubic  cardinal  spline  with  function  values  equal  to 


zero  apart  from  the  zero  point,  where  it  assumes  the  value  L,  is  given 
by 

«0 

S(x)  =  /3  (-2  +  B(x-j).  (2.4- 

Its  behavior  can  be  judged  from  the  graph  shown  in  Figure  2.4: 


X 

Fig.  2.4  Fundamental  Cardinal  Cubic  Spline 


This  spline  which  is  known  as  "fundamental  cardinal  cubic  spline"  has 
unbounded  support,  is  twice  continuously  differentiable,  consists  of 
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cubic  polynomials  within  each  interval,  is  symmetric  and  interpolates 

the  infinite  data  vector  [...,0,0,  1,0,0,...].  Each  data  vector 

different  from  the  above  one  can  uniquely  be  spline -interpolated  by  a 

linear  combination  of  fundamental  splines  with  the  function  values  at  the 

grid  points  as  coefficients, 

00 
«— 1 

F(x)  =  i  fk  S(x-k)  . 

In  the  discussion  above  we  have  limited  ourselves  to  cardinal  cubic 
spline  interpolation  on  the  whole  real  line.  This  is  for  many  reasons 
a  restriction: 

a)  in  all  practical  problems  data  are  given  on  a  limited 
part  of  the  real  line; 

b)  often  the  data  are  not  regularly  distributed  with  grid 
distance  =  1; 

c)  cubic  spline  interpolation  is  just  one  kind  of  spline 
interpolation;  why  not  use  another  one,  say  quadratic 
or  quintic  ? 

There  is  a  good  reason  to  discuss  very  isolated  cases:  because 
they  show  up  the  very  behaviour  of  the  gbneral  solution.  This  is  why 
cardinal  spline  interpolation  has  been  discussed  and  not  spline  inter¬ 
polation  on  an  irregular  limited  grid  of  data.  In  such  cases  the  formulas 
are  nc  longer  as  simple  (the  matrix  in  (2.4-6)'  is  not  as  regular,  but 
still  tri-diagonal)  --  the  features,  however,  remain  the  same.  As 

far  as  point  (c)  is  concerned,  there  is  a  simple  answer:  cubic  spline 
interpolation  has  been  chosen  not  only  as  a  compromise  between  linear 
interpolation  and  spline  interpolation  of  highest  smoothness  (which  is  a 
sinx/x  -  interpolation  as  shown  in  (H.  Slinkel,  1980)  );  it  has  also 

been  chosen  for  serious  practical  reasons:  the  cubic  spline  still 
retaines  a  high  degree  of  simplicity,  while  its  attractive  features 
(smoothness,  localness)>  usually  adherent  to  more  sophisticated  interpola- 


tion  functions,  remain.  Two  of  its  properties  are  worth  being  at  least 
mentioned:  the  cubic  spline  minimizes  the  overall  squared  second 
derivatives  (which,  in  the  case  of  small  first  derivatives,  approximately 
corresponds  to  a  minimization  of  the  overall  curvature  and,  therefore, 
the  elastic  energy), 

J  [F"(x)l2  dx  =  min.  (2.4-9) 

among  all  possible  interpolating  functions.  This  property  is  called  the 
minimum  norm  property.  The  second  one,  called  best  approximation 
property,  guarantees  that  the  interpolating  spline  has  smaller  distance 
from  a  given  function  (sampled  at  the  data  points)  than  any  other  non¬ 
interpolating  spline;  distance,  in  this  context  is  defined  via  the  pseudo 
norm  (2.4-9). 

As  far  as  approximation  properties  are  concerned,  the  following 
error  estimates  can  be  shown  to  hold  (Ahlberg  et  al.,  1967): 

|F^(x)-f^(x)|  *  h”®  J|D4f(x)|dx,  a  =  0,  1  (2.4-10) 

with  h  denoting  the  grid  distance  and  £(x)  the  function  to  be  approximated 
by  the  cubic  spline  F(x) .  ,  Similar  error  bounds  hold  for  second  and  third 
derivatives. 


2.4.2  Two-dimensional  cubic  spline. 


Analogous  to  the  one-dimensional  case  one  can  define  a  function 
of  the  independent  variables  x  and  y,  interpolating  all  data  on  an  infinite 
regular  rectangular  grid  with  constant  grid  distance  equal  to  1  (cardinal 
grid)  such  that  the  interpolating  function  is  twice  continuously  differ¬ 
entiable  with  respect  to  both  independent  variables  x  and  y: 


F(x,y) 


c  kl  B(x-k)  B(y-l). 


(2.4-11) 


* 

( 
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As  in  the  one-dimensional  case,  the  coefficients  {  c  ^  ]  are  not  simply 
the  function  values  at  the  grid  points,  but  linear  combinations  of  function 
values  with  derivatives.  The  grid  information  is  hardly  ever 
available  and,  therefore,  one  intends  to  find  base  functions  such  that  the 
function  values  at  the  grid  points  are  identical  with  the  coefficients 
fc  in  (2.4-11).  This  is  in  fact  possible;  the  unique  fundamental 
bicubic  cardinal  spline  function  is  just  a  product  of  the  one -dimensional 
splines  discussed  in  Chapter  2.4.  1: 


S(x,  y)  = 


w  mo 

l  l 


L  a  .a.  B  (x-i)  B  (x-j). 

-CO  1  J 


(2.4-12) 


This  spline  which  is  known  as  "fundamental  cardinal  bicubic  spline"  has 
unbounded  support,  is  twice  continuously  differentiable  with  respect  to 
x  and  y,  consists  of  bicubic  polynomials  within  each  cell  of  the  grid, 
and  interpolates  the  infinite  data  array  6^  6^  : 

i  •  » 

•  •  • 

«  •  • 

...0  0  0... 

...0  1  0... 

...  0  0  0. . . 

•  •  • 

•  •  • 

•  •  • 

Its  behavior  is  similar  to  that  of  the  one-dimensional  spline  . 

Figure  2.5  shows  its  main  features. 

Each  data  array  which  is  different  from  the  above  one  can  uniquely  be 
spline-interpolated  by  a  linear  combination  of  such  fundamental  splines 
with  the  function  values  at  the  grid  points  as  coefficients, 


F(x,  y) 


w  w 

.11 


L  L  ft1  S(x-k)  s(y-l). 

«  kl 


(2.4-13) 


As  before  we  have  limited  ourselves  to  the  presentation  of  cardinal 
cubic  splines  defined  on  the  whole  two-dimensional  plane.  This  ideal 
case  will  never  be  met  in  practical  applications;  a  bounded  support 


Fig.  2.5  Fundamental  bicubic  cardinal  spline 

causes  the  formulas  to  become  slightly  more  complicated,  but  the  main 
features  are  retained.  A  non-uniform  grid  is  also  possible,  however, 
it  has  to  be  generated  by  lines  parallel  to  the  coordinate  lines.  GSPP 
does  not  deal  with  such  a  case. 

As  far  as  minimum  properties  are  concerned,  there  is  also 
an  analogy  with  the  one-dimensional  spline)  the  bicubic  spline 
minimizes  the  following  integral: 


(2.4-14) 


a4  F(x,y)  • 
3x2  S  y2 


2 

dxdy  =  min. 


among  all  possible  interpolating  functions  (minimum  norm  property) 
and  the  interpolating  spline  has  least  distance  (where  "distance"  is 
defined  via  the  norm  (2. 4- 14))  from  a  given  function,  sampled  at  the 
grid  points,  than  any  other  non-interpolating  spline  (Ahlberg  et.  al. ,  1967). 

As  far  as  approximation  properties  are  concerned,  it  can  be  shown 
that  the  approximation  error  depends  on  the  grid  spacing  in  the  following 
way: 


dyF(x,  y)  _  Syf(x,  y) 
dx*  bx?  dx“  3xe 


3-a  3-a 

+hy  ), 


y  =  a+fl  s  6  ,  ^  3,  jis  3 

and  ,  h^  grid  distances  in  x  and  y-direction. 

2. 5  The  frequency  content  of  a  bicubic  spline  surface 

A  bicubic  spline  surface  is  based  on  data  (function  values) 
distributed  on  a  regular  rectangular  grid.  Formulated  differently,  the 
bicubic  interpolating  spline  is  a  smooth  function  interpolating  all  samples 
of  the  original  function. 

It  is  well  known  from  the  sampling  theorem  (see  e.g. 

E.O.  Brigham,  1974,  p.  83  ff)  that  the  original  function  can  be 

exactly  reconstructed  from  the  samples,  only  if  the  original  function  is 

baud -limited  with  highest  frequency  ^max»  fche  sampling  interval  h 

is  smaller  than  or  equal  to  1/2  f  ;  the  frequency  1/h  =  2f  is 

n  max  max 

called  the  Nyquist  sampling  rate. 

In  general,  the  sampled  functions  are  not  band-limited  and  no 
interpolation  function  is  able  to  exactly  reproduce  the  original;  this  fact 
is  called  aliasing. 
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To  know  which  frequencies  can  be  represented  by  a  cubic  inter¬ 
polating  spline  is  of  interest  by  itself;  moreover,  since  operations  like 
differentiations  are  performed  not  on  the  data  but  on  the  interpolated  values, 
it  is  important  to  know  how  the  interpolated  values  respond  under  such 
operations;  this  can  be  seen  best  in  the  frequency  domain. 

In  order  to  study  the  frequency  behavior  of  bicubic  splines,  we 
remember  that,  according  to  equation  (2.4-13),  the  two-dimensional 
spline  base  function  is  just  a  product  of  one -dimensional  splines. 

Therefore,  we  investigate  now  the  frequency  behavior  of  the  one¬ 
dimensional  spline. 

2.5.  1  Spectrum  of  the  cubic  spline.  We  define  the  spectrum 

of  a  function  f(x)  as  its  Fourier  transform 

•0 

F(w)  =  J £(x)e  ^  dx  ,  i  =  VTf  (2. 5- la) 

with  its  inverse 

•0 

f(x)  =  ~  jF(o))eit0Xdu)  .  (2. 5- lb) 

We  are  interested  in  the  case  of  f{x)  to  'be  a  cubic  spline.  There  are 
at  least  two  ways  of  approach:  a  delicate4  and  mathematically  rather 
involved  one  which  starts  with  the  spectrum  of  the  B- spline,  takes  ad¬ 
vantage  of  the  properties  of  Euler- Frobenius  polynomials,  and  derives 
the  transform  of  fundamental  cardinal  cubic  splines;  this  approach  car 
be  found  in  detail  in  (H.  SUnkel,  1980).  Here  a  much  simpler  and 
straight  forward  derivation  will  be  given  which  is  based  on  formulas 
derived  in  (H.  Sunkel,  1977a).  In  the  following  it  will  be  assumed  that 
data  are  given  at  all  integers  from  to  +®  and  that  a  cubic  spline 

should  be  fitted  to  these  data.  Furthermore,  the  data  should  be  such  that 
the  integral  (2.5  -la)  exists.  Since  we  know  from  Chapter  2.4  that  the 
cubic  spline  consists  of  cubic  polynomials  within  each  interval  (between 
two  consecutive  data),  we  can  split  up  the  integral  (2.  5- la)  in  a  countable 
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infinite  number  of  integrals,  each  one  being  extended  over  one  interval 
of  length  1.  Consequently,  the  contribution  of  the  cubic  polynomial 
defined  between  the  integers  m  and  m+1  is  given  by 
m+1 

(m),  >  _  f  (m),  -icox 
F  '(<o)  =  J  f  (x)  e  dx 


m+1 

P  Am),  .  .  ,(m),  . 

=  I  f'  '(x)  cos  cox  -  x  f  '(x) 


sin  cox 


(2.5-2) 


(in) 

with  the  cubic  polynomial  f'  (x).  Taking  (2.  4-. . . )  into  account, 
equation  (2.5-2)  can  be  split  up  further  into 


F(m(l)  .  L  4m)  •  is<m)U] 

k=0 


(2.  5-3a) 


with  being  polynomial  coefficients,  and  cj^^co)  and  S^m\co) 

the  real  and  imaginary  part  of  the  k  ^  degree  polynomial  contribution 
for  the  interval  [m,  m+l]: 


(2.  5-3b) 


the  result  of  these  simple  integrations  can  be  found  in  (H.  Siinkel,  1977a, 
p.  37  ff. ): 


c(m) 

m+1 

cos  to  x 

k 

J  (x-m) 

c  (m) 

m 

sin  co  x 

. 

„(m),  .  .Jm),  .  1  r  -io)(m+l)  -icon. 

Co  \Co)  -  iSo  '(co)  =  -r-[e  -e  ] 


j  -ico  (m+i) 

t—  e  , 

ioo 


k  =  1,2,... 


(2.5-4) 
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After  some  manipulations  it  is  possible  to  express  the  contribution  of 
a  cubic  polynomial,  between  the  interval  [m,  m+i  ],  to  the  spectrum 
F(co)  in  the  form 


F(m)Q  =  t  aH'*"1" 

I _ 


- kTT  tl  (2-5-5) 

(ia)  k+‘  r=0  (k'r)- 


The  four  polynomial  coefficients  (a,^11  3  ,  k”0,...,3,  however,  depend 

linearly  on  the  function  values  f  ,  f  ,  and  the  derivatives  f'  ,  f' 

m  m+i  m  m+j 

The  goal  is'  ;:o  figure  out  all  coefficients  of  f  ,  f  ,  f  ,  and  f1 

m  m+i  m  m+* 

and  finally  to  express  f1  and  f1  ,  by  the  infinite  vector  {f  }, 

7  r  m  m+i  7  mJ  ' 

m  =  , . . , 

Evaluating  (2.5-5)  explicitly  ,  one  finds  the  "contribution"  of 


(2.  5 -6  a) 


all  f  to  be 
m 

•0 

7  f  e“iwm 

r  12  / .  -li 

f  (1-e 

V-l 

L  m 

'  60 

-•D 

00 

r*i 

•0 

since  Li  .. 

m+i 

-icom 

ico  y 

e  — 

e  L 

the  contribution  of  all  f  can  be  shown  to  be 

m+i 


«D 


to 


(2.5-6b) 


In  a  similar  way  one  obtains  the  contributions  of  f'  and  f'  ,  : 

'  m  m+i 

l  V'iamt  (2.5-5C) 


and  1  f'  4  4  (Z+e^  )  . 


W  £0 


(2.  5-6d) 


Introducing  ecu+e  =  2  cosco  and  e  w  -  e  =  2i  •  sin  w  ,  the 
sum  of  (2.5-6a)  +  (2.5-6b)  and  (2.5-6c)  +  (2.5-6d)  can  be  written  as 


w 

V  -item  r  24  ,  .12  .  , 

Li  e  *  l  774  (1“cos£o)“773  sm  &)  ], 


m 


Ol 


09 


(2.  5-7 a) 


i  • )  ,.  -icom  r12  4  .  Nl 

lj  f'  e  L — 5  sm  CO  -  —  (2+cos  U!)J. 


m 


0) 


or 


(2.  5-7b) 


There  remains  still  £'  to  be  expressed  by  f  : 

m  1  m 


it  can  be 


shown  that  £'  is  a  linear  combination  of  f. 

m  J 


f'  =  -3  I  a.  (f  ..  -  f  .) 
m  j  =  i  J  m+J  m“3 


(2.5-8) 


with  coefficients  Qi.  =  a5  =(-2+/3)^.  Taking  this  relation  into  account 

3 


7.  fi  e  lwrn  can  be  transformed  into 
m 


E,,  -icon  „  )  -icom  )  /,  „  . 

f'  e  =  -3 L  e  u  a.  (f  ,.-f  .)  ; 

“  j  =  1  1  m+j 


m 


m=-« 


taking  advantage  of  the  identity 

CO 

I 


m--“ 


eluma.£ 

J  3 


m=-'° 


m 


which  is  the  frequency  equivalent  of  a  "time- shift",  and  interchanging 
the  sequences  of  summation  over  m  and  j,  one  obtains 


50  w? 

If  e-ium  =  [-3Za.(eij“e-ii“)Hl 


.  -ioxn-i 
f  e  J  , 


and  with  -  e  =  2i  •  sin  (j'to), 

•o  m  oo 

Z  f'  e  =  -6i[  Za.  sin(jto)]*[  Z 

„  m  j  _ 

3  =  1  m= 


-itom, 
f  e  1. 

m  J 


(2.5-9) 


Because  |  a|  =  0.2768  <  1,  the  first  sum  in  above  equation  can 
be  expressed  in  closed  form  (Gradshteyn,  No.  1.447,  p.  40), 


V  4  V  a  sin  to 

L  o’  sin(ju)  =  L 

.1  =  1  J=1 


which,  with  a  =  -2  +  / 3  reduces  to 


•0 

I  mow)"-]  • 


(2.5-10) 


and,  with  (2.5-7b),  (2.5-9),  and  (2.5-10),  the  contribution  of  all  f^ 
to  the  spectrum  is  given  by 


,  36  sin2t0  ,  12  \  Y  *  -iwm 

( - + - sin  (t) )  •  L>  f  e 

to4(2+ cos  to)  to3  m 


Adding  the  contribution  of  f^  from  equation  (2.  5 -7 a),  we  obtain  for  the 
whole  spectrum 


,  12  (1-cos  to)2  Y  ,  -itom 

F(to)  =  — - \  L  f  e 

to  (2+  cos  to  )  _«  m 


(2.5-11) 


with  the  identity 
4  (l-cosfc))2 


.  to 

sm-  4 


(  -nr1-  ) 
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the  Fourier  transform  of  the  cubic  spline  is  finally  given  by 
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sin 
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F(a>)  = 


2  +  cos 


a) 
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1 


L  e-iWm 

„  m 


(2.5-12) 


This  result  is  very  interesting  and  deserves  to  be  discussed: 
Let  us  first  investigate  the  summation  part  of  F(co): 


Fs(to)  =  I 


f  e 
m 


•i60m 


(2.5-13) 


This  is  precisely  the  discrete  Fourier  transform  of  the  infinite  data 

vector  (f  }  ,  m  =  -*>,  ...,®  (see  e.g.,  E.O.  Brigham,  1974).  The 
m 

fore,  its  inverse  Fourier  transform  should  again  be  the  data  vector: 


&  ■  s  J 


00  00 

£  1  £  J  «'* 
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(2.5-14) 


The  integral  in  the  above  formula  can  also  be  written  as 

CO 

m  p 

r  -ico  m  icox  ,  ioo(x-m)  , 

J  e  e  dco  =  J  e  dw 


_  00 


=  [cos  6o(x_m)dai  +  i  J  sinco(x-m)  dco  . 


The  second  integral  vanishes  because  of  the  asymmetry  of  the  integrand. 
The  first  integral  can  be  split  up  into 
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we  conclude  from  the  convolution  theorem  that  a  cardinal  cubic  spline 
f(x)  corresponds  to  F (to)  with  f(x)  being  the  result  of  a  convolution  of 
a  fundamental  cardinal  cubic  spline  S(x)  with  the  data  sequence 
(sequence  of  impulses), 


f  (x)  =  S(x)  *  *  L  i  6  (x-m) 

-CO 


CO  #0 

.L  I 


S(t)  6  (x-m-t)dt. 


f  S(x-m) 
m 


The  above  equation  is  identical  with  (2.4-4)  and  so  is  the  circle  closed. 

Now  we  can  also  seethe  more  elegant  way  of  deriving  the  Fourier  trans¬ 
form  of  the  cubic  spline:  find  the  Fourier  transform  of  the  fundamental 
cardinal  cubic  spline  and  multiply  it  with  the  discrete  Fourier  trans¬ 
form  of  the  data  sequence.  This  approach  can  be  found  in  (H.  SUnkel,  1980). 

Let  us  once  more  consider  F^  (<*')  and  assume  that  F^, (to)  =  0  for 
|co!>  IT  ,  or  in  other  words,  let  F  (to)  be  the  Fourier  transform  of  a 

U 

frequency  band  limited  process  (function).  Then  the  sum  in  equation 
(2.5-15)  reduces  to  a  single  element,  the  integral 


IT  .  , 

f  ,  \ ,  ,  sm  IT  (x-m) 

cos  to  (x-m)dto  =  2  - 1 

x-m 


which,  with  (2.5-14),  leads  to  the  interpolation  function 


f(x)  =  Zj  f 


sin  TT(x-m) 


m  IT  (x-m) 


(2.5-16) 


This  function  is  interesting  insofar  as  it  can  be  shown  to  be  an  inter¬ 
polating  spline  function  of  highest  possible  degree"  (I.J.  Schoenberg,  1973). 
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It  is  infinitely  often  continuously  differentiable,  but  suffers  from  localness. 
It  is  also  remarkable  in  that  it  is  exactly  the  function  which  reproduces 
an  original  IT-  band -limited  function  sampled  at  a  sampling  rate  equal 
to  1.  This  is  the  essence  of  the  sampling  theorem.  Therefore,  using 
(2.5-16)  as  interpolation  function,  corresponds  to  the  assumption  of  an 
originally  band-limited  function  with  highest  frequency  co  =  TT  (the 
Nyquist  frequency). 

The  cubic  spline,  however,  shows  a  somehow  different  behavior: 
strictly  speaking,  its  spectrum  is  not  band-limited,  buc  the  Fourier 
transform  of  the  fundamental  cardinal  cubic  spline 


3 

2  +  cos  co 


sin 


co 


co 

1 


■)4 


(2.5-17) 


is  such  that  it  practically  annihilates  all  frequencies  above  co  =  27T  . 
This  factor  is  of  interest  because  it  varies  with  the  degree  of  the 
spline.  In  the  case  of  a  step  function,  which  is  the  spline  of  lowest 
degree  --0,  this  dampening  factor  degenerates  to 


S0(co) 


(2.5-18) 


in  the  case  of  the  function  (2.5-16),  which  is  the  spline  of  highest 
degree  --  *,  this  factor  degenerates  to  a  window 


S«(co) 


1  for  |  co  |  ^  V 
0  otherwi.'e 


(2.5-19) 


In  between  these  two  extrema  there  is  the  large  family  of  splines  of  all 
possible  degrees.  A  graph  of  the  factor  (2.5-17),  which  is  at  the 
same  time  the  Fourier  transform  of  the  fundamental  cardinal  cubic 
spline,  is  shown  in  Fig.  2.6. 
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frequencies  u)  >  it  .  The  reason  for  this  is  that  the  integral  of  the 
Fourier  transform  of  the  cardinal  spline  (and  ail  other  splines  as 
proven  in  (H.  Sunkel,  1980)  )  is  equal  to  IT, 


go 

,-s 


?,+  COS  O) 


.  40 

sin  I 

40 
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)4  d  40 


ff. 


Consequently,  the  amount  of  dampening  in  the  low  frequencies  40  *  17 
is  equal  to  the  build-up  of  frequencies  for  40>  TT  : 
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2  +  cos  40 
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Summarizing  we  can  say,  that  the  frequency  dampening  of  a  spline 
of  arbitrary  degree  in  the  frequency  range  jtoj  ^  TT  is  exactly  com¬ 
pensated  in  the  range  (  TT,  «°)  ,  relative  to  the  window  (2.5-19)*  The 
higher  the  degree  of  the  spline,  the  better  is  the  approximation  to  this 
window.  Therefore,  the  dampening  of  the  lower  frequencies  and  the 
build-up  of  higher  frequencies  must  be  caused  by  the  limited  degree 
of  continuous  differentiability  of  the  spline  (a  kind  of  Gibb's  phenomenon). 
For  this  reason  we  conclude  that  the  frequencies  40 >  TT  are  not 
reliable  anymore  --  it  is  more  or  less  frequency  noise  which  is  only 
present  in  order  to  compensate  the  deficiency  in  the  range  40 &  IT  . 

To  understand  this  is  essential  in  order  to  get  a  better  idea  of 
how  reliable  differentiated  splines  are. 

J 

2.5.  1.1  Spectrum  of  the  differentiated  cubic  spline.  The 
module  GSPP  is  also  capable  of  providing  derivates  of  profiles  and 
surfaces  up  to  and  including  second  order.  In  order  to  better  under¬ 
stand  the  reliability  of  the  output,  which  is  a  differentiated  cubic 
(bicubic)  spline,  it  is  essential  to  investigate  the  impact  of  a  differentia¬ 
tion  on  the  spline  in  the  frequency  domain. 
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This  time  we  go  the  short  way  and  make  use  of  the  fact  that  the 
cardinal  cubic  spline  is  the  result  of  a  convolution  of  a  fundamental 
cardinal  cubic  spline  with  an  impulse  comb  of  data  (the  sequence 

tt  J ), 

wi  • 


£(x)  =  S(x)  *  L  f  6  (x-m)  # 

_ m 


(2.5-20) 


The  differentiated  cardinal  cubic  spline  is  likewise  given  by 


f'(x)  =  S'(x) 


fm  6  (x-m) 


(2.5-21) 


The  corresponding  equivalent  in  the  frequency  domain  is  a  product  of 
the  Fourier  transform  of  S'(x)  with  the  Fourier  transform  of  the 
impulse  comb;  the  latter  is  given  by  equation  (2.5-13), 
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-  ito  m 


From  (2.5-12)  we  know  the  Fourier  transform  of  S(x)  to  be 


.  to 

sm  — 

- - -  ( - —  )4 

2  +  co  s  to  '  to 


and  therefore,  S(x)  can  be  obtained  from  the  inverse  Fourier  transform, 


S (x)  =  —  J  S(to)eiaJX  dto 


S'(x)  as  derivative  of  S(x)  with  respect  to  x  is  then  given  by 


s'(x)  =  J  S(to)  to  e1W  X  dto  , 


(2.5-22) 


and  its  Fourier  transform  is  simply 
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s(>)(55  =  ±  ]  J  s j 

since  we  know  from  equation  (2.5-15)  that 

CO  _ 

f  ±i(o>  -  co  )x  ,  __  e/  — . 

Je  '  '  dx  =  2TT  •  6  (o)  -  0)  )  , 

-  •o 

we  obtain 

00 

S^(co)  =  i  J  S{  w)  o)6(o)-o)  )do)  » 

„CO 

which,  because  of  the  property  of  the  6-  function,  reduces  to 

S^(co)  =  i  co5 (to)  .  (2.5-23) 

And  finally,  the  Fourier  transform  of  the  differentiated  cardinal  cubic 
spline  is 

F(i)(co)  =  S(J)(co)  •  Fs(to)  , 
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sm  -  4 

( - -) 

0) 

2 


l 


.  -iO)m 
f  e 
m 


(2.5-24) 


We  recognize  the  well-known  fact  that  a  differentiation  dampens  the 


lower  (co  <  1)  and  amplifies  the  higher  (o) >  1)  frequencies.  For 
comparison  purposes  we  also  give  the  Fourier  transform  of  the  highest 
possible  degree  ( 50 )  cardinal  spline, 


03 

iw_% 
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for 
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• 

The  graphs  in  Fig.  2.7  give  a  comparison  of  both  transforms  for  first 
and  .  ocond  order  derivatives.  Notice  the  shift  of  the  frequency  sensitivity 
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FOURIER  TRANSFORM  OF  THE  1.  DERIVATIVE  OF 

FUNDAMENTAL  CARDINAL  SPLINES  OF  ODD  OEGREE 

•K  ...  TRANSFORM  OF  THE  1.  DERIVATIVE  OF  SPLINE  OF  OEGREE  K 


U 

FOURIER  TRANSFORM  OF  THE  2.  DERIVATIVE  OF 

FUNDAMENTAL  CARDINAL  SPLINES  OF  GOD  OEGREE 

«K  ...  TRANSFORM  OF  THE  1.  DERIVATIVE  OF  SPLINE  Or  OEGREE  K 


u 

Fig.  2.7  Fourier  transforms  of  spline  derivatives 
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maximum  towards  u i  =  it  with  increasing  number  of  differentiations. 
What  are  the  consequences  of  this?  Since  the  frequency  resolution 
becomes  worse  for  higher  frequencies,  the  reliability  of  spline 
differentiation  decreases  with  each  differentiation  and  the  frequency 
noise  represented  by  the  frequency  part  co  >  IT  gains  more  influence 
on  the  result  of  the  differentiation. 

2.5.2  Spectrum  of  the  bicubic  spline.  Recalling  the  defining 
equation  (2.4-13)  of  the  bicubic  spline 


00 

[,y)=  I  I 


f  S(x-m)  S(y-n)  , 

mn  '  3  ' 


(2.5-25) 


with  [f  }  the  infinite  data  array  and  S(x)  the  fundamental  cardinal 
cubic  spline,  it  is  obvious  that  the  bicubic  spline  is  a  tensor  product 
of  cubic  splines  in  x  and  y.  Therefore  it  is  pretty  easy  to  find  the  two- 
dimensional  Fourier  transform 


F(a)x,W  )  =  J  J  f(x,y)e'l(C°xX  +  4°yy)  dxdy  . 


(2.5-26) 


Since  (2.5-25)  is  a  two-dimensional  convolution  of  a  two-dimensional 
impulse  comb  with  a  two-dimensional  fundamental  cardinal  cubic  spline 


f(x,  y)  =  S(x)  S(y)  *  * 


n 


f  6  (x-m)  6(y-n), 
mn 


(2.5-27) 


the  corresponding  Fourier  transform  F (to  x,  w  y)  is,  according  to  the 
convolution  theorem,  simply  a  product  of  the  Fourier  transform  of  the 
2-D  fundamental  cardinal  cubic  spline  and  the  Fourier  transform  of 
the  2-D  impulse  comb. 

To  find  the  first  is  easy,  it  is  just  a  product  of  two  one- 
dimensional  transforms, 

S  (wx,  W  )  =  3(&»x)  S(coy). 
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S(cu  ,  to  )  = 
x  y 


.  ojv  .  y 

_ 9 _ ,Sm2  <Sin  2  4 

(2+cos  COx)  (2  +  cos  (jiy)  0)x  oy  ’ 


(2.5-28) 


the  latter  is  simply  the  2-D  analogue  to  (2.5-13), 


FS(tOx,  Wy) 


=  L  1  f  +  “V1) 


(2.5-29) 


.And  finally,  with  (2.5-27),  the  2-D  Fourier  transformation  of  the 
cardinal  bicubic  spline  is  given  by 
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(2.5-30) 


Expressed  in  terms  of  the  bicubic  spline  coefficients,  the  above  equation 
has  the  form  which  is  analogous  to  (2.5-5)  of 

“  “  1  1 
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(2.5-31) 
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An  equivalent  expression  can  also  be  found  in  (Bhattacharyya,  1969)* 

All  that  has  been  stated  so  far  for  the  one -dimensional  spline,  carries 
over  to  the  Uvo -dimensional  one: 

There  exists  a  highest  possible  degree  (  ",  ®)  2-D  cardinal  spline, 
which  is  fully  analogous  to  (2.5-16),  with  the  form 


/  \  ^  sin  IT  (x-m)  .  sin  If  (y-n) 

(x,  y)  =  L  L,  £  — - r—  — L-r—4  . 

nm  ff  (x-m)  ft  (y-n) 


(2.5-32) 


This  function  is  infinitely  often  differentiable,  but  suffers  from  localne-s 
if  compared  with  the  bicubic  spline.  It  is  the  interpolating  function 
which  reproduces  a  band-limited  function  with  highest  frequencies 
(0X  =  If ,  u)y  =  if  at  a  sampling  rate  equal  to  1.  Therefore,  using 
(2.5-32)  as  interpolation  function,  corresponds  to  the  assumption  of  an 
originally  band-limited  function;  if  this  was  not  the  case  (as  usual  in 
almost  all  practical  applications),  the  interpolating  function  suffers  from 
aliasing  effects. 

The  bicubic  spline's  spectrum  is,  strictly  speaking,  not  band- 
limited,  but  the  2-D  Fourier  transform  of  the  fundamental  cardinal 
bicubic  spline  (2. 5-28) anniliates  practically  all  frequencies  above 
=  2TT,  ojy  «  2ir.  The  frequencies  above  oj  =  ff  ,  however,  are 
caused  by  the  discontinunity  of  the  splines  third  derivative.  It  is  a 
kind  of  Gibb's  phenomenon;  therefore,  these  higher  frequencies  represent 
more  or  less  only  "discontinuity  noise",  which  is  exactly  compensated 
in  the  frequency  range  |u)  |  s  if  .  This  compensation  causes  the  lower 
frequencies  to  be  dampened  (compare  Fig.  2.7  ). 

2.6  Practical  1-D  and  2-D  spline  routines 

2.6.1  The  cubic  spline  routine.  In  chapters  2.4  and  2.5  a 
somewhat  artificial  case  has  been  discussed:  splines  based  on  an 
infin  te  number  of  data,  uniformly  distributed  on  a  grid  with  constant 
grid  distance  equal  to  1.  The  discussion  has  been  performed  for  this 


I 


exotic  case  because  it  is  considerably  simpler  to  point  out  the  essential 
features;  also  the  finite  spline  behaves  likewise  because  of  its  localness. 
For  applications,  however,  it  is  necessary  to  deal  with  the  general 


case. 

Detailed  derivations  can  be  found  in  a  standard  text  on  splines 

like  (Ahlberg  et  al.,  1967).  Here,  only  the  very  necessary  things  will 
* 

be  presented. 

Let  there  be  given  a  sequence  of  abscissas  {xm}  ,  m=l  , 

...,  M;  f  =  f(xm)*  An  interpolating  cubic  spline  is  a  twice  continuously 
differentiable  function  over  the  range  [a,  b]  which  consists  of  (M-l)  cubic 
polynomials,  each  of  them  defined  on  an  interval  [xj,  x-+J  ]  , 


(*) 


(2.6-1) 


The  coefficients  {a^  }  ,  m=l  .  .  .  ,  M-l,  can  be  found  via  the  continuity 

conditions  for  first  or  second  order  derivatives.  This  leads  to  a  system 
of  M  linear  equations 

Af"  =  Bf  (2. 6 -2a) 


A  is  a  strongly  diagonally  dominant  tridiagonal  symmetric  matrix  with 
positive  diagonal  elements;  it  follows  from  well  known  theorems  in 
matrix  algebra,  that  A  is  positive  definite.  The  uniqueness  of  the  spline 
is  guaranteed  if,  in  addition  to  the  function  values,  the  second  derivatives 
are  given  at  the  end  points  a  and  b  of  the  interval  [a,  b].  Since  such 
kind  of  boundary  informations  are  hardly  ever  available,  the  most 
natural  choice  is  to  assume  them  to  be  zero.  This  corresponds  to  a 
linear  behavior  of  the  spline  at  the  end  points  a  and  b  of  the  interval. 

The  corresponding  cubic  spline  is  referred  to  as  a  "natural  spline". 

GSPP  assumes  vanishing  2nd  derivatives  at  the  boundary  points,  and 
therefore,  calculates  exclusively  natural  cubic  splines. 
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The  solution  of  (2.6-2) 
f"  =  A-*Bf. 


(2. 6  -2b) 


can  be  found  very  efficiently  by  a  Gauss-Jordan  elimination  procedure 
(Spath,  1973;  p.  10  ff. ).  If  the  system  has  been  solved,  it  is  fairly 
easy  to  find  the  coefficients  Ca^m^3,  k  =  0,  3;  m=l,  ...,  M- 1: 
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(f'L +  2f"  ) 


(2.6-3) 


(f «  .  -  f «  ) 
m+i  m  ' 


with  Ax  :  =  x  ,  -x  and  Af  :  =  f  -  f 
m  m+i  m  m  m+l  m 

In  a  similar  way  (2. 6 -2a,  b)  and  (2.6-3)  can  be  expressed  in 
terms  of  1st  derivatives,  which  will  be  used  in  the  calculation  of 
bicubic  spline  coefficients. 

2.6.2  The  bicubic  spline  routine.  Let  a  regular  rectangular 

grid  consist  of  M«N  gridpoints  and  let  the  function  values  at  the  grid 

points  be  {f  3,  m=l,  ...,  M;  n=l,  ...,  N.  Then  a  bicubic  spline 
r  mn  r 

consists  of  (M-1)(N-1)  bicubic  polynomials  in  x  and  y 


3  3 


(m,n),  .  Y  )  (m,  n),  >k,  «l 

E  '(x,  y)  L  L  a,.  '(x-xm)  (y-y) 


k=0  1=0 


m'  w  yn 


(2.6-4) 


ftn  n) 

with  coefficients  ’  )  expressed  by  the  product 

A(m,n)  =  t>Kn)j  .  ) 


(2.6  -5a) 


***m!J  +~*&*r**-  ■W-i'"  -  V 


and  the  grid  spacing  matrices  H(h^)  and  H(h  ), 


1  0  -3/h2 

0  1  -2/h 

0  0  3/h2 

0  0  -  1/h 


2/h3 

1/h2 

■2/h3 

1/h2  j 


(2.6 -5b) 


where  h  =  grid  distance  in  x-direction 
x 

h^  =  grid  distance  in  y-direction 


F  is  a  matrix  with  data  information 


F  = 


f 

mn 

q  f  , 

mn  m,  n+) 

^m,  n+l 

p 

mn 

r  p 

mn  rm,  n+l 

rm,  n+i 

Xm+»,  n 

£ 

qm+i,n  m+j,  n+l 

^m+i,  n+l 

p 

m+i,  n 

rm+i,n  ^m+i ,  n+l 

rm+i,  n+i 

(2. 6-5c) 


where 

f  =  function  value  at  the  gridpoint  (m,n) 
mn  ° 

p  =  first  x-derivative  at  the  gridpoint  (m,  n) 
mn 

q  =  first  y-derivative  at  the  gridpoint  (m,  n) 
mn 

r  =  second  xy-derivative  at  the  gridpoint  (m,  n). 
mn 

The  derivatives  {p  }  ,  {q  )  ,  and  [r  }  at  all  gridpoints  are 

mn  mn  mn 

determined  by  continuity  conditions  of  second  order  derivatives  similar 

to  the  one-dimensional  case.  A  unique  bicubic  spline  representation, 

however  requires,  apart  from  all  function  values  [f  }  at  the  grid 

mn 

points,  the  following  additional  boundary  informations: 

5-Pj,  n^  «  fejvl  ’  n“^’  •••» 

’  (2.6 -5d) 

.  Um>N)  »  m=l  ,  ....  M; 
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{r  3  ,  m=l,  M;  n=l,  N. 
mn 

Since  these  boundary  informations  (derivatives)  are  hardly  ever 
available  in  practical  applications,  similarly  as  in  the  one- dimensional 
case,  assumptions  have  to  be  made  concerning  their  values.  The  most 
natural  assumption  is  to  assume  vanishing  second  order  derivatives  along 
the  boundary  normal.  This  is  an  implicite  assumption  relative  to  the 
boundary  data  (2.6-5d).  The  practical  determination  of  the  coefficients 
runs  then  as  follows: 


1. )  Solve  the  spline  equations  (2.6-2)  for  first  order  deriva¬ 

tives  as  unknowns  with  data  vectors  {f  3,  m=l, ....  M 

mn 

and  vanishing  2nd  derivatives  at  the  boundary  points, 

D2f(x,  y  )|  _  =D2f(x,  y  )|  =0  ,  for  all  "columns" 

x  n  X"X j  x  n  x=xj^j 

n=l . N  obtaining  so  all  1.  derivatives  in  x-direction 

[p _ 3,  m=l,  ....  M;  n=l,  ...,  N. 

mn 

2. )  Interchange  the  role  of  x  and  y:  Solve  the  spline  equations 

(2.6-2)  for  first  order  derivatives  as  unknowns  with  data 

vectors  {f  3,  n=l,  ...,  N  and  varnishing  2nd  derivatives 
mn 

at  the  boundary  points,  D2f(x  ,y)|  =D2f(x  A  =0, 

y  '  m  y=yj  y  m,y)'y=yN 

for  all  "rows"  m=l,  ...,  M  obtaining  so  all  first  derivatives 

in  y-direction  {q  3,  m=l,  ...,  M;  n=l,  ...,  N. 
mn 

3. )  Replace  "f"  by  "q"  and  solve  according  point  (1)  or,  which 

is  equivalent,  replace  "f"  by  "r"  and  solve  according  point 
(2).  The  result  is  edl  second  order  mixed  derivatives 

1  r  3  ,  m—  1 ,  . . . ,  j  n-1,  ...,  ISf. 

mn 

With  these  values  obtained,  the  bicubic  spline  is  completely  defined. 

(m  n) 

(The  coefficients  (a,  *  3  depend  only  on  {f  ,  p  ,  q  ,  r  3, 

Ki  mn  mn  mn  mn 

m=l,  ...,  M;  n=l,  ...,  N  and  on  the  grid  spacing  h  and  h  .) 

x  y 
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Therefore,  it  is  necessary  to  have  this  set  of  defining  values,  the 

quadrupel  (f,  p,q,  r)  at  each  grid  point,  available  in  order  to  define 

mn 

a  spline  uniquely.  Generally  speaking,  also  the  grid  distances  h  and 

h  can  change  with  m,  h  =  h  (m),  h  can  change  with  n,  h  =  h  (n). 
y  x  x y  y  yv 

However,  GSPP  always  assumes  h^  and  h^  to  be  constant  (they  need 
not  necessarily  be  equal).  This  simplifies  computations  quite  consider¬ 
ably:  the  grid  information  matrices  become  constant;  therefore,  it  is 
possible  to  perform  all  calculation  normalized,  with  a  constant  grid 
spacing  in  x  and  y-direction  equal  to  1.  Of  course,  the  calculated 
derivatives  have  to  be  interpreted  accordingly  --  they  refer  to  grid 
spacing  equal  to  1  and  have  to  be  scaled  by  h^  and  h  later  on  in  order 
to  obtain  real  values.  The  gain  is  two-fold:  firstly,  all  derivative 
calculations  (solutions  of  the  spline  equations)  are  highly  stabilized,  and 
secondly,  the  calculation  of  the  product  (2.6-5a)  --  which,  in  GSPP,  is 
not  formulated  in  terms  of  two  matrix  products  but  is  programmed 
explicitly  --  gains  in  calculation  speed  because  a  number  of  divisions 
(or  multiplications) is  avoided  ,  This  sounds  trivial  but  is  very  essential 
when  the  grid  is  extended  and  if  a  huge  number  of  interpolations  has 
to  be  performed. 

2.6.3  Interpolation /differentiation  of  splines.  Interpolation 
with  cubic  and  bicubic  splines  is  really  trivial  as  soon  as  the  defining 
values  are  available  --  function  values  and  1st,  or  2nd  order  deriva¬ 
tives  at  the  knots  {x  ),  m~l,  ...,  M  in  the  one-dimensional  case, 

m 

and  the  quadruple  {f,  p,q,  r)  ,  m=l,  ...,  M;  n=l,  ...,  N  in  the 

mn 

two-dimensional  case.  Therefore,  an  interpolation  of  a  function  value 
at  a  point  x  (x,  y)  requires  the  following  steps  to  be  performed: 

1)  find  the  interval  number  m  (grid  element  number  s  m,  n)  to  which 
the  point  x  (x,  y)  belongs; 

2)  calculate  the  cubic  spline  coefficients  using  (2.6-3)  (bicubic 
spline  coefficients  using  (2.6-5)  ; 


) 
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3)  perform  the  product  (2.6-1),  for  the  cubic,  (2.6-4)  fc.r 
the  bicubic  spline. 

Computation  efficiency  can  be  gained  by  performing  the  products 
(2.6-1)  and  (2.6-4)  not  blindly,  but  by  reducing  the  number  of  operations 
involved  to  a  minimum  like 


f(x)  =  a0  +  x  (a,  +  x  (a2  +  x  a3)  )  ,  x  =  x-x^  (2. 6 -6a) 

for  the  cubic  spline,  and 

f(x)  =  aoo  +  y  (aoi  +  y  (a  02  *  y  a03)  ) 

+  x  ((al0  +  y  (au  +  y  (aJ2  +  y  al3))) 

(2.6-7a) 

+  x  ((a20  +  y  (a2l  +  y  (a22  +  y  a23))) 

+  x  (a30  +  y  (a3j  +  y  (a32  +  y  a33)))))  ,  x  =  x-Xm>  y  =  y-yn 


for  the  bicubic  spline.  (This  arrangement  reduces  the  number  of  multi¬ 
plications  from  originally  at  least  28  to  15  and  keeps  the  number  of 
additions  constant.)  If  the  coefficients  refer  to  grid  spacing  1  (normal¬ 
ized  coefficients),  then  x  ^x-x^J/h^  and  y  =  (y-y^/h^  . 

Any  derivative  0s  (X  ^  3  of  the  cubic  spline  (2.6-1)  can  be 
expressed  by 


D^f(x) 

x 


l  - 
\  <  <*> x 


k-a 


0  s  a  £  3; 


k-a 


(2.6-8) 


however,  it  is  more  efficient  to  explicitly  write  down  the  derivatives: 
f'(x)  =  a-i  +  x  (2a2  +x.  3a3), 

f"(x)  =  2a2  +  6a3x,  (2.6-6b) 

f'"(x)  =  6a3 

Similarly,  any  derivative  0  =»  0^  s  3,  0  s  a2  s  3  of  the  bicubic  spline 

(2.6-4)  can  be  expressed  by 
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oc  ,  5 

D  f(x,y)  = 


« 


Sx^y*2 


=  «,!«,!  1  l  'k''1-k'a-1'°!! 


-  -  aki  C )  C  ^  y 

k=ttj  i=ofe  w  tt2 


(2.6-9) 


with  the  double-index  a:  =  (ott,  a2),  |«|  =  fll,+  a,,  0  s  fl,,  a2  s  3. 

This  compressed  expression,  however,  is  far  from  optimal  as  far  as 
CPU-time  saving  is  concerned.  Therefore,  all  partial  derivatives  should 
rather  be  used  in  the  following  forms  (see  also  (Sunkel,  1980)  ): 

fx(x,y)  =  aio  +  y  (au  +  y  (aB  +  y  ai3  ))  (2.  6-7b) 

+  x  (2(azo  +  y(a2J  +  y  (a22  +  y  a23))) 

+  x  •  3  (a30  +  y  (a3t  +  y  (a32  +  y  a33)))), 

fy.(^,y)  “  a0|  x  (all  x  (a21  x  a3j)) 

+  y  (2(a02  +  X  (a, 2  +  x  (a22  +  x  a32))) 

+  y*3  (a03  +  x  (al3  +x  (a23  +  x  a33)))), 

£xx(x>  y)  =  2(a20  +  y  (a2l  +  y  (a22  +  y  a23)) 

+  x  •  3(a30  +  y  (a31  +  y  (32  +  y  a33)))), 

fXy(x»  y)  aU  *  y  (2aJ2  +  y  3a13) 

+  x  (2(a21  +  y  (2a22  +  y*3a23)) 

+  x*  3  (a3J  +  y  (2a32  +  y.  3a33))), 

fyy(x»  y)  =  2(a02  +  x  (a, 2  +  x  (a22  +  x  a32  )) 

+  y  •  3  (a03  +  x  (al3  +  x  (a23  +  x  a33)))), 

fJDDI(x>y)  =  6(a3o  +  y  (a3j  +  y  (a32  +  y  a33)))» 
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f 

xxy 


(*.  y) 


f 

xyy 


(x,  y) 


2(a2l  +  y  £a22  +  y  •  3  a23) 

+  x  •  3(a31  +  y  (2a32  +  y*3a33))), 
2(aJ2  +  x  (2a22  +  x-  3a32) 


+  y*  3  (al3  +  x  (2a23  +  x*3a33))), 


f  (x,  y)  =  6(a03  +x  (at3  +  x  (a23  +  x  a33  ))), 

yyy 

*„r(x»y)  =  6(a3»  +  y  (2a32  +  y*3a33)), 

f  (x,  y)  =  4(a22  +  y  3a23  +  x*  3  (a32  +  y3a33)), 

xx  yy 

f  (x,  y)  =  6(%3  +x  (2a23  +  x‘3a33))  , 

xyyy 

f  (x,  y)  =  12  (a32  +  y*  3a33), 

xxxyy  1 

f  (x,  y)  =  12  (a23  +  x*  3a33), 

xxyyy 

f _  (x,  y)  =  36  a33  . 

xxxyyy 


GSPP  is  capable  of  providing  plots  of  derivatives  0  ^  Oq  ,  a2  5  2. 

2.7  C  onto  ur  finding 

Let  the  prediction  of  function  values  at  a  regular  rectangular 
grid  be  done  and  assume  that  the  spline  routines  have  generated  the 
quadruple  (f,p,  q,  r]^,  m=l,  , M;  n=l,  . ..,  N  of  bicubic  spline 
defining  values.  In  other  words,  let  a  bicubic  spline  surface  be 
given  such  that  on  each  grid  element  Jxm*  xm+1  »  yn»  the 

smooth  interpolation  function  given  by  (2.6-4).  The  goal  is  to  find 

the  contour  f(x,  y)  =  constant  =  c. 

Consider  just  one  single  grid  element  R  with  the  bicubic  poly- 
J  mn 

nomial  (2.6-4)  as  interpolating  element.  Then  it  is  obvious  that  at 
least  one  contour  exists  if  the  condition 

min  f(x,  y)  <  c  <  max  f(x,  y) 

(x'y)fRmn  (x-y  >eRtan 

is  fulfilled.  In  the  case  of  one-sided  equality  the  contour  degenerates 
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to  a  point,  in  the  case  of  two-sided  equality  the  bicubic  element  is 
flat  and  horizontal,  all  coefficients  {a^}  »  k,l  =  0,  ...,  3  are  zero 
apart  from  a00  (which,  naturally,  can  also  vanish),  andno  contour 
exists. 


The  decision  process  of  existence  or  non-existence  of  a  contour 
(or  more  contours)  is  only  the  beginning  of  the  long  procedure  of  actually 
finding  them.  Remember  that  we  have  to  find  the  solution(s)  of  the 
equation 


ffr,  y) 


I  k  1 

L  aklX  y 


k=0  1=0 


=  c  =  const. , 


(2.7-1) 


a  slightly  simplified  form  of  (2.6-4).  In  one  dimension  this  corresponds 
to  the  solution  of 

f(x)  =  a<j  +  ajx  +  a2x2  +  a3x3  =  c, 
which  is  nothing  else  but  finding  the  zeros  (  s  3)  of 

(a<)-c)  +  a3x  +  a2x2  +  a3x3  =  0  (2.7-2) 

with  almost  arbitrary  coefficients.  Even  in  such  a  simple  case,  a 
simple  solution  usually  does  not  exist.  The  approximation  methods 
used  for  the  solution  of  above  cubic  equation  essentially  consist  of 
finding  the  intersection  of  the  horizontal  line  fj(x)  =  c  and  a  line 
segment  of  (2.7-2)  which  represents  the  behavior  of  this  cubic  poly¬ 
nomial  in  the  neighborhood  of  the  zero  being  considered.  In  other  words, 
the  zeros  can  be  found  by  approximating  the  cubic  polynomial  (2.7-2)  by 
a  continuous  and  piecewise  linear  function,  and  intersecting  this  approxi¬ 
mation  function  with  the  line  f3(x)  =  c.  It  is  obvious  that  each  line  seg¬ 
ment  of  the  piecewise  linear  function  has  either  1  intersection  or  none 
(Fig.  2.7.  1) 

In  view  of  these  facts  it  should  be  clear  that  direct  contour¬ 
finding  with  the  bicubic  function  (2.7-1)  is  hopeless.  There  is,  however, 
a  way  to  overcome  this  problem,  similar  to  that  discussed  above: 
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Fig.  2.7.1  Intersection  of  f(x)  with  (x) 


By  approximating  the  bicubic  function  (2.7-1)  by  a  number  of  simple 
elements  which  allow  for  a  simple  contour -finding. 

The  simplest  such  approximation  function  is  obviously  a  plane. 

A  plane  either  has  1  intersection  with  another  (horizontal)  plane  or  has 
no  intersection;  moreover,  the  intersection  is  a  straight  line.  A  plane 
is  uniquely  defined  by  3  parameters  which  can  be  taken  as  the  function 
values  at  3  not  coinciding  points.  Therefore,  one  would  conclude  that 
piecewise  linear  (flat)  triangular  elements  are  the  obvious  choice  for  the 
approximation  of  any  continuous  two-dimensional  function  analogous  to 
the  one-dimensional  case. 

A  linear  triangular  element  is  defined  by  the  equation 


f(x,  y)  =  b0  +  bjx  +  b2y 
and,  intersected  with  the  horizontal  plane 


(2.7-3) 


(2.7-4) 


gives  an  intersection 

y(x)  =~  (c  -  b0  -  bjx) 

which  is  obviously  the  equation  of  a  straight  line.  Since  the  element  is 
flat,  only  two  boundary  lines  of  the  triangular  element  can  have  an 
intersection  with  (2.7-4)  or  no  boundary  line  has  an  intersection.  To 
find  these  intersections  is  indeed  very  simple:  first  it  is  necessary  to 
check  if 

min  f.  <  c  <  max  f-  , 

i— 1,2,3  i=  1,2,3 

where  f .  ,  i  =  1,2,3  are  the  function  values  at  the  3  corners  of  the 
triangle;  the  second  step  consists  in  the  actual  calculation  of  the  2 
intersection  points.  The  straight  line  connecting  these  two  intersection 
points  is  part  of  a  contour.  The  full  contour  is  then  the  continuous 
and  piecewise  linear  line  consisting  of  the  above  described  line  segment 
(Fig.  2.7.2). 


Fig.  2.7.2  Contours  on  piecewise  linear  triangular  elements 
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This  is  in  principle  the  contour-finding  procedure  frequently 
used  in  connection  with  detailed  digital  terrain  models:  sampled  terrain 
heights  and  other  kinds  of  structure  informations  like  rivers,  lakes, 
roads,  dams,  artificial  structures  are  taken  into  account,  a  dense 
"triangulation  is  generated",  and  the  terrain  heights  at  the  corner  points 
of  the  triangles  are  either  known  or  are  predicted.  It  should  be 
pointed  out  that  all  this  sounds  very  simple,  but  is  extremely  complicated. 
The  enormous  problem  consists  primarily  of  analyzing  how  the  human 
intellect  "produces"  an  image  on  the  basis  of  data  and  secondly,  to 
"translate"  this  stream  of  logical  operations  in  a  programming  language. 

We  are  not  going  to  discuss  this  method  in  detail  because  GSPP 
has  been  designed  for  another  possible  method  of  contour  finding. 

Recall  that  the  bicubic  polynomial  (2.7-1),  as  part  of  the  bicubic 
spline  surface,  is  defined  on  a  rectangle.  Therefore  it  is  quite  natural 
to  approximate  the  surface  not  by  flat  triangular  elements,  but  by  simple 
elements  also  defined  on  a  rectangle.  The  simplest  element  defined  on  a 
rectangle  is  a  function  with  4  parameters  which  can  be  determined  from 
the  function  values  at  the  4  corners  of  the  rectangle.  Such  a  function, 
in  general,  is  no  plane  anymore;  it  is  a  hyperbolic  paraboloid  (saddle 
surface)  with  the  equation 

f(x,  y)  =  b0  +  bjx  +  b2y  +  b3xy  (2. 7 -5a) 

this  function  is  bilinear,  which  can  be  seen  more  easily  in  the  equivalent 
form 

f(x,  y)  =  (c0  +  Cjx)(d0  +  dty)  (2. 7 -5b) 


CodO 

~  b0, 

C|d0 

~  ul> 

Codi 

=  bz. 

ctdj 

n 

43 

II 

or  in  the  form 
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1 

V 

f(x,  y)  =  L 
k=0 

with 

boo  =  bo» 
bto  =  bi* 
b0»  =  b2> 
b„  =  b3. 

Bilinear  means  that  the  function  is  linear  along  the  coordinate  lines 
x  =  const,  and  along  the  coordinate  lines  y  =  const.  In  any  other 
direction  the  function,  in  general,  not  linear  (Fig.  2.7.3). 

This  bilinear  element  plays  a  central  role  in  contouring  based 
on  rectangles  and  a  thorough  discussion  is  essential  for  an  understand¬ 
ing  of  the  whole  contouring  logics.  We  will,  therefore,  investigate  its 
properties  in  detail. 

Let  us  first  introduce  a  new  coordinate  system  (x',y')  parallel 
to  the  old  system  (x,  y)  with  origin  coordinates  (x0,  y0)  in  order  to 
eliminate  the  linear  terms  in  (2.7-5a), 

x  =  x'  +  x0,  y  =  y'  +  y0  •  (2.7-6) 

Then  the  intersection  between  the  bilinear  element  (2.7-5a)  and  a  hori¬ 
zontal  plane  fj(x,  y)  =  c  =  constant  assumes  the  form 

c  =  b0  +  bt(x‘  r  x0)  +  b2(y'  +  y0)  +  vb3(x'  +  x0)(y'  +  y0) 

=  (b0  +  bjx0  +  b2y0  +  b3x0y0)  +  (bj  +  b3y0)x' 

+  (b2  +  b3x0)y'  +  b3x'y'  ;  (2.7-7) 

the  linear  terms  (x1,  y')  vanish  if  b3  +  b3y0  =  0  and  b2  +  b3x0  =  0; 
these  conditions  provide  the  coordinates  of  the  origin  of  the  new 
coordinate  sy stem  (x1,  y1) 


1 

N7  k  1 

L  b  x  y  (2. 7 -5c) 

1=0  Ki 


x0 


.  y o  = 


(2.7-8) 


Figure  2.7.3  :  A  bilinear  element 


With  these  values  equation  (2.7-7)  takes  on  the  simple  form 


x'y'  = 


(c-b0  + 


(2 


Since  the  right  hand  side  of  the  equality  sign  is  constant,  the  line  of 
intersection  has  the  equation  x'y1  =  constant.  But  this  is  exactly  the 
equation  of  a  hyperbola  with  asymptotic  lines  parallel  to  the  coordinate 
lines  (because  x'  =  const. /y'  and  y'  =  const. /x').  Therefore,  we  know 
also  the  axes  of  the  hyperbola:  they  are  mutually  orthogonal  and  span 
an  angle  of  45°  with  the  coordinate  lines.  In  order  to  prove  this  we 
introduce  a  new  coordinate  system  (x,  y)  rotated  by  an  angle  a  : 


=  x  cos  ft  -  y  sin  ft 


rl  = 


x  sin  ft  +  y  cos  ft 


(2.7-10) 


With  this  coordinate  transformation,  the  product  x'y'  assumes  the  form 


x'y'  =  (x2  -ya  )  cos  ft  sin  ft  +  x  y(cos2  ft  -  sin2  ft  ) 


=  (x2  -y 2  )—  sin  2ft  +  xy  cos  2  ft  =  const. 


The  above  equation  becomes  purely  quadratic  if  cos  2  ft  =  0  which  corres- 
7T 


ponds  to  a  =  —  proving  the  statement  made  above.  So  we  finally  obtain 


the  equation  of  the  contour 


x2  -  y2 


r*  (c  -  bo  +  ) 


(2.7-11) 


This  purely  quadratic  expression  is  the  mid-point  equation  of  a  hyperbola. 
It  refers  to  a  coordinate  system  (x,y)  which  is  shifted  relative  to  the 
original  system  (x,  y)  by  (x0,  y0)  having  values  (2.7-8),  and  rotated  by  an 
angle  of  45°.  The  axes  of  the  hyperbola  coincide  with  the  new  coordinate 
axis  (x  =  0,  y  =  0)  and  span  therefore  the  same  angle  of  45°  with  the 
original  system.  The  hyperbolas  as)  .nptotes  are  parallel  to  the 
original  coordinate  lines  (Fig.  2.7.4). 

After  having  pointed  out  all  relevant  facts  concerning  the  line  of 
intersection,  we  come  back  now  to  the  bilinear  element.  The  bilinear 
element  is  linear  along  each  coordinate  line;  therefore,  it  is  also  linear 
along  the  4  boundary  lines  of  the  rectangle.  Each  linear  function  can 
be  intersected  at  most  once  by  another  linear  function,  and,  for  the 
same  reason,  can  there  be  at  most  one  contour  intersection  point  cn 


y=y0  (asymptote) 


Figure  2.7.4  Contours  of  a  bilinear  element 


each  of  the  4  boundary  lines  of  the  rectangle.  Before  we  discuss  this 
maximum  case  we  state  that  the  other  two  possibilities  are  2  intersection 
points  or  no  intersection  at  all.  In  the  case  of  only  two  intersection 
points  it  is  clear  that  both  points  belong  to  one  and  the  same  hyperbola 
(because  the  asymptotic  lines  are  parallel  to  the  original  coordinate 
lines),  and  the  connection  between  the  two  points  is  clear.  The  way 
of  making  the  connection  between  intersection  points  is  not  so  evident 
(Schumaker,  1976,  p.  249  ff).  Is  there  a  simple  and  unique  answer  to  how 
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the  4  intersection  points  in  Fig.  2.7.5  have  to  be  connected? 


P2 

Figure  2.7.5:  Bilinear  elements  --  contour /boundary- 
intersection  points  . 


Yes,  there  is  a  simple  answer  (  Sunkel,  1977a).  The  clue  is 
that  the  asymptotic  lines  of  the  hyperbolas  are  parallel  to  the  coordinate 
lines.  Therefore,  the  position  of  1  out  of  the  4  intersection  points 
relative  to  the  origin  (x0,  y0)  of  the  new  coordinate  system  determines  the 
way  of  point  connection  uniquely:  In  Fig.  2.7.5  the  point  Pj  is 
below  the  asymptotic  line  x  =  x0;  therefore,  it  can  only  be  connected 
with  Pz  whose  coordinate  y  is  smaller  than  y0  ;  since  the  hyperbolas 
are  symmetric  relative  to  the  origin,  the  connection  of  P3  with  P4  is 
automatically  fixed.  Consequently,  the  simple  calculations  of  the  origin 
(2.7-8)  provide  sufficient  information  about  the  way  of  point  connection 
(Fig.  2.7.6). 


Figure  2.7.6:  Bilinear  elements  --  connection  of  contour/ 
boundary  intersection  points. 


The  coordinates  of  the  hyperbolas  origin  are  valuable  also  for 
another  reason:  they  tell  if  4  intersection  points  are  possible  (2  hyper¬ 
bolas)  or  if  just  2  intersection  points  (1  hyperbola)  exists  within  the 
square;  if  the  origin  is  located  outside  the  square,  then  at  most  one 
hyperbola  exists,  if  it  is  located  within  the  square,  two  hyperbolas  are 
possible.  There  is  an  exceptional  case  which  should  also  be  mentioned: 
the  case  of  hyperbolas  degenerating  into  the  asymptotic  lines;  this 
happens  if  the  x-coordinates  of  one  pair  of  points  (and  the  y-coordinates 
of  the  other  pair  of  points)  coincide  --  such  lines  are  known  as  saddle 
lines  (see  Fig.  2.7.6). 

All  the  considerations  made  above  about  the  intersection  curve(s) 
and  the  way  of  intersection-point  connection  are  simple  but  essential 
for  the  logics  of  contour  finding. 

Let  us  go  back  to  the  coefficients  {b^L  k  =  0,  ...,•  3,  in 
equation  (2. 7 -5a).  They  can  easily  be  determined  as  linear  combinations 


of  the  function  values  at  the  4  corners  uf  the  square: 


b0  =  foo 

bj  =  fjo  "  ^oo 

^2  ~  foi  “  *oo 

^3  =  tfjl  ■  ^to)  -  (foi  "  ^oo) 


(2.7-12) 


with  f^  =  f(x  =  k,  y  =  1).  Here  it  has  been  assumed  that  the  square 
on  which  the  particular  bilinear  element  is  defined,  has  unit  length. 

This  saves  a  number  of  calculations  (divisions  or  multiplications)  in 
GSPP. 

Summarizing,  the  process  of  contour -finding  for  a  bicubic  spline 
surface  consists  basically  of  the  following  steps: 

1)  Interpolate  the  bicubic  spline  at  the  grid  points  of  a 
rectangular  array  which  is  a  proper  subdivision  of 
the  original  grid  on  which  the  spline  has  been  defined 
(Fig.  2.7.7).  The  function  approximating  the  bicubic 
spline  is  defined  by  the  array  of  function  values  at  these 
subgrid  points.  The  function  is  continuous  and  con¬ 
sists  of  bilinear  elements,  each  defined  on  one  individual 
subgrid  element.  Store  this  approximating  function. 

2)  Start  with  the  lowest  possible  contour  value  and  search 
for  a  '-ontour  intersection  point  along  the  boundary 
rectangle  of  the  whole  subgrid  --  contours  intersect¬ 
ing  the  boundary  are  open  (non-periodic)  with  respect 
to  the  contouring  area  in  consideration.  Having  found 
the  first  intersection  points,  (here  are  two  alternatives 
to  find  the  next  one  in  this  element:  either  by  first 
comparing  the  contour  value  with  the  other  4  remaining 
contour  values,  singling  out  the  possible  line(s)  and 
calculating  the  coordinates  of  the  intersection  point(s), 
or  to  take  advantage  of  the  relation  between  the  hyper¬ 
bolas  and  their  asymptotic  lines.  Since  the  origin  coordinates 


Figure  2.7.7:  Grid  (4x5)  and  subgrid  (16x21) 


are  necessary  anyway  for  the.  decision  of  how  many  hyperbolas 
are  possible  and  for  the  point  connection,  a  hybrid  solu¬ 
tion  has  been  chosen.  This  intersection  point  searching 
continues  like  a  domino  from  one  element  to  the  next 
until  the  contour  leaves  one  of  the  four  boundaries. 

A  certain  integer  array  associated  with  the  array  of 
bilinear  elements  "remembers"  the  position  of  all  inter¬ 
section  points  relative  to  the  boundary  rectangle  of  a 
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particular  element,  and  "directs"  the  program  part 
which  is  responsible  for  the  calculation  of  the  contour  inter¬ 
section  points.  This  is  necessary  in  order  to  avoid  a 
multiple  calculation  (and  plotting)  of  contours. 

After  all  non-periodic  contours  (for  a  particular 
conto'ur  value)  have  been  found,  the  program  searches 
for  periodic  contours  (contours  which  do  not  cross  the 
boundary). 

3)  Plotting  of  the  contour.  Repeat  steps  (2)  and  (3)  until 
the  contour  value  exceeds  the  maximum  function  value  of 
the  surface. 

These  are  the  basic  steps  which  are  necessary  to  find  contours 
fo  a  bicubic  spline  surface.  However,  there  are  and/or  can  be  a 
number  of  secondary  procedures  involved  which  make  the  steps  (2)  and 
(3)  rather  complicated.  In  the  sequel,  three  intermediate  procedures 
are  described  which  GSPP  is  capable  of  performing. 

2.8  Optional  contour  procedures 

We  start  here  with  the  most  commonly  used  procedure,  the 
contour  smoothing. 

2.8.1  Contour  smoothing.  Theoretically,  a  contour  smoothing 
would  not  be  necessary  if  the  hyperbolas  of  the  last  section  (contours 
of  a  bilinear  element)  could  be  sampled  at  a  sufficiently  high  rate. 

This,  however,  is  relatively  expensive  in  terms  of  CPU-time.  There¬ 
fore,  it  has  been  decided  to  take  only  2  points  of  a  single  hyperbola 
into  account,  its  2  intersection  points  with  the  boundary  of  the  square  on 
which  the  bilinear  element  is  defined.  The  actual  (non-normalized) 
size  of  the  square  (subgrid  distance)  is  chosen  as  approximately  2mm 
by  GSPP  and  can  be  changed  within  certain  bounds  by  the  user.  How¬ 
ever,  consecutive  contour  points  will  have  a  mutual  distance  of  up  to 
approximately  fl*  subgrid  distance  (»  3.5  mm  in  this  case).  A 
linear  connection  of  these  contour  points  might  give  a  too  rough  picture. 
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Therefore,  it  is  either  necessary  to  decrease  the  subgrid  distance 
or  to  smooth  the  contours.  GSPP  will  always  smooth  the  contours 
unless  the  smoothing  is  suppressed  by  the  user.  The  interpolation 
function  used  for  the  smooth  representation  of  the  contours  is  a 
partitioned  and  overlapping  cubic  parameter  spline. 

A  cubic  parameter  spline  is  a  cubic  spline  whose  coordinates 
x  and  y  depend  on  a  parameter  which  is  here  taken  to  be  an  approxi¬ 
mation  of  the  arc  lengths  of  the  contour, 

x  =  x(s) 

y  =  y(s). 

Therefore,  both  x  =  x(s)  and  y  =  y(s)  are  cubic  splines  depending  on  the 
parameter  s.  In  GSPP,  the  arc  length  is  taken  as  the  accumulated 
length  of  the  piecewise  linear  function  determined  by  the  originally 
calculated  contour  points. 

The  number  of  points  for  a  single  contour  may  become  quite 
large  (a  couple  of  hundred);  moreover,  the  points  are  by  no  means 
equally  spaced.  These  circumstances  can  cause  instabilities  in  the 
spline  algorithms.  In  order  to  avoid  this  it  has  been  decided  to  first 
"clean-up"  the  crude  vector  of  contour  points  which  is  to  be  understood 
as  an  elimination  of  all  contour  points  whose  mutual  distance  is 
smaller  than  some  percentage  of  the  subgrid  distance  (assumed  to  be 
25%  but  can  be  changed  by  the  user  within  some  bounds).  After  this 
clean-up  process,  the  vector  of  contour  points  is  subdivided  into  a 
number  of  overlapping  vectors  (GSPP:  25  points/vector,  6  points  overlapping). 
The  overlapping  has  been  chosen  in  order  to  preserve  the  smooth  transi¬ 
tion  from  one  part  of  the  contour  to  the  next  part  to  such  a  degree  that 
it  is  not  possible  to  detect  the  transition  by  visual  means;  a  6-point 
overlapping  fulfills  this  requirement;  the  number  25  has  been  chosen 
for  reasons  of  stability. 

After  the  calculation  of  the  spline  parameters,  further  contour 
points  are  interpolated  such  that  the  overall  mutual  distance  of  the 
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contour  points  varies  only  between  narrow  bounds.  Fig.  7 
comparison  between  smoothed  and  non- smoothed  contours. 


Figure  2.8.  1 


rn  i  t  Tvnn 77  ^v\\\\wv  v  ^ 


.8.  1  gives  a 


Smoothed  and 
non-  smoothed 
contours 
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2.8.2  Contour  mapping.  Usually,  the  two-dimensional  coordin¬ 
ates  x  and  y  are  interpreted  as  cartesian  coordinates  and  the  mapping 
"data  space  “*  image  space"  (contour  plot)  is  assumed  to  be  given  by 

x*  =  c  •  x 

*  -  .  (2.8-1) 

where  (x,  y)  are  the  coordinates  of  a  data  point,  (x*  ,  y*)  are  the 
coordinates  of  the  plotted  point,  and  c  is  the  constant  scale  factor. 

However,  especially  in  all  earth  sciences  the  coordinates  are 
not  to  be  understood  as  cartesian  coordinates  but  as,  e.g.,  spherical 
coordinates  9,  A,  and  there  is  the  wish  or  need  to  choose  a  mapping 
different  from  that  in  equation  (2.8-1).  For  example,  a  Mercator 
projection 

x*  =  C|  In  [tan  (^  +  )]•  sign  (9) 

y*  =  c2  (A-A0) 

with 

9,  X  =  spherical  coordinates  (latitude,  longitude), 

A0  =  longitude  origin 

C|,  c2  =  scale  factors  in  x  and  y-  direction; 

Another  example,  a  linear  transformation  of  the  form 

x#  =  C|  (x  cos  a  +  y  sin  Q() 
y*  =  c2  (-x  sin  a  +  y  cos  a) 

with  the  constant  azimuth  a  as  rotation  angle.  Obviously,  there  is  an 
unlimited  variety  of  possibilities  for  coordinate  transformations 
(mappings).  If  any  mapping  different  from  (2.8-1)  is  desired,  the 
user  has  to  define  this  mapping  by  providing  the  corresponding  sub¬ 
routine  and  he  has  to  inform  GSPP  that  a  mapping  is  requested.  This 
mapping  has,  naturally,  to  be  done  for  all  points  of  the  contour; 
actually,  it  is  performed  after  the  clean-up  of  the  contour  point 
vector.  Fig.  2.8.2  shows  a  Mercator  projection  of  the  contours  in 
Fig.  2.8.  1. 


Fig.  2.8.2  Mercator  projection  of  the  contours  shown  in  Fig.  2.8.1 


2.8.3  Contour  labeling.  By  labeling  of  contours  we  under¬ 
stand  the  plot  of  the  contour  value  into  an  interval  which  is  kept 
free  of  the  contour  line.  The  contour  has  to  have  a  sufficient  length 
such  that  a  contour  value  fits  in  an  interval  which  is  smaller  than 
the  actual  length  of  the  contour.  If  a  contour  is  too  short,  no  label 
(contour  value)  will  appear.  The  direction  of  the  label  is  designed  to 
be  identical  with  the  tangent  to  the  contour  at  the  midpoint  of  the  label 
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The  length  of  the  label  depends  on  the  actual  contour  value,  the 
number  of  decimal  places  to  be  plotted,  and  the  sign  of  the  contour 
value,  according  to  the  equation 

(  n+2+u(-x)  for  |x|  ^  10.0 


w  =  c*h* 


(2.8-2) 


n+2  +  int(logl0{  |x  j ))  +  u(-x)  for  |x|  >  10.0 


w  =  label  length, 

h  =  symbol  height, 

c  =  ratio  symbol  length/ symbol  height 
x  =  contour  value, 

n  =  number  of  decimal  places  to  appear, 

int  =  integer  function 

u(x)  =  unit  step  function,  u(x)  = 


1  for  x>  0 
0  for  x  s  0  * 


If  not  defined  explicitly  by  the  user,  the  number  of  decimal 
places  is  chosen  such  that  at  least  3  significant  digits  appear.  If 
the  absolute  contour  value  is  greater  than  or  equal  to  1000,  no  decimal 
digit  and  no  decimal  point  will  appear.  If  the  contour  value  is  zero, 

A 

only  a  "0"  without  decimal  point  and  decimal  digit  will  be  plotted  (see 
Fig.  2.8.  1) . 


Moreover,  labeled  contours  can  be  plotted  with  multiple  line- 
width  (if  the  plotter  is  designed  for  such  a  purpose). 

The  labeled  contour  interval  has  to  be  an  integer  multiple  of  the 
non-labeled  contour  interval  ;  the  default  value  in  GSPP  is  each 
second  contour  labeled  with  double  linewidth. 

2.8.4  Contours  within  a  window.  Usually,  the  complete  bicubic 
spline  surface  from  m=l,  ...,  M  and  n=l,  ...,  N  will  be  contoured. 
Sometimes,  however,  one  is  probably  interested  in  only  a  part  of  it 
in  order  to  see  more  details  when  plotted  in  a  larger  scale  or  one 
likes  to  have  more  contours  than  usual  plotted  in  a  particular  area. 

For  such  purposes  GSPP  offers  the  possibility  to  plot  only  a  rectangular 
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This  feature,  however,  is  limited  in  many  respects:  it  allows 
on’y  rectangular  windows  on  full  grid  elements  (from  one  grid  point  to 
another  grid  po*  <t).  Furthermore,  it  is  not  possible  to  suppress  the 
contour  drawing  within  a  rectangular  region.  In  order  to  overcome 
all  these  deficiencies,  a  highly  sophisticated  routine  has  been  developed 
and  incorporated  into  GSPP;  it  will  be  described  in  the  following  section. 

2.8.5  Contour-free  regions  of  arbitrary  shape.  Sometim e s 
there  is  the  desire  to  plot  the  contours  only  within  a  predefined  closed 
region,  or  vice  versa,  to  suppress  the  contour  drawing  within  certain 
regions.  For  example,  it  might  be  requested  to  plot  a  gravity  anomaly 
..lap  only  for  the  state  of  Ohio  or,  a  rather  recent  application, 
to  generate  a  world  geoid  solely  based  on  altimeter  data  and  to  plot 
this  geoid  only  over  the  oceans  and  greater  lakes  because  of  the  lack  of 
data  over  the  continents  (e.g.,  Rapp,  1979)*  Before  we  describe 
the  procedure  which  is  capable  of  doing  this,  let  us  mention  that  such  a 
statement  is  very  simple  to  make  --  its  transir.^on  into  a  computer 
language,  however,  is  extraordinarily  complicate*.  The  reason  is 
primarily  that  such  a  region  is,  in  general,  not  conve..  practically 
arbitrary  in  shape,  and  last  but  by  no  means  least,  may  consist  of  a 
very  huge  number  of  boundary  points.  In  order  to  give  an  idea:  the 
detailed  world  shore-outline  data  bank  consists  of  some  80  000  points. 
Moreover,  there  might  be  contour-regions  within  contour-free  regions, 
and  so  on  --  this  could  be  the  case  in  the  above  mentioned  altimeter 
example:  geoid  plot  only  over  "water":  ocean  (continent  (lake 
(island  (lake  •))))••  More  technically,  the  procedure  should 

be  able  to  hariul  sequences  of  subsets 

A)"3  Aj  3  A2  3  . . .  3  Afl  (2.8-3) 


such  that  contours  are  drawn  only  on 
Asfe  “  Aa^-j  k  =  0,  1,  ...,“ 


(2. 8-4a) 
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pr  only  on  its  complement 

A?^- j  —  Agj^  >  k  —  1»  2,  ~  (2. 8-4b) 

(see  Fig.  2.8.4  ).  Sometimes  it  is  even  required  that  the  contours 
are  drawn  (or  suppressed)  on  a  union  of  mutually  exclusive  (disjoint) 
sequences  of  subsets.  In  this  case  A^  and  n  depend  on  a  parameter  i 
(see  Fig.  2.8.4): 

Aq  Aj  . . .  — *  A*  ,  i  *  1,  . . . ,  I.  (2. 8-5 ) 

ni 

i 

GSPP  is  capable  of  handling  even  this  case. 

Honestly  speaking,  the  author  himself  had  absolutely  no  reasonable 
idea  how  to  start  when  he  got  involved  in  this  problem.  It  is  clear  that, 
in  order  to  plot  contours  only  in  certain  regions,  it  is  necessary  to 
a)  find  the  intersection  between  the  contours,  and  b)  to  know  which 
part  of  the  contour  has  to  be  plotted  and  which  one  has  to  be  deleted. 

In  order  to  find  the  intersections  between  a  contour  and  all  region 
boundaries  it  is  necessary  to  develop  a  stable  algorithm  which  is  capable 
of  finding  the  intersection  between  two  straight  lines,  and  to  run  this 
algorithm  for  all  possible  line  segments  of  the  contour  and  all  possible 
line  segments  of  the  boundaries.  If  the  number  of  boundary  points  is 
very  small,  this  procedure  will  not  be  very  expensive;  however,  if  the 
boundaries  consist  of  100  points  and  more,  this  method  becomes 
absolutely  prohibitive  in  terms  of  computer  time.  How  else  should 
one  attack  this  problem?  Since  in  nature  almost  everything  is 
optimized  (very  probably  because  of  the  long  time  of  evolution),  it  came 
to  the  author's  mind  to  analyze  the  almost  automatic  human  decision 
process  for  the  particular  problem  of  suppressing  the  drawing  of  contours 
within  predefined  regions,  and  to  "translate"  this  process  step  by  step 
in  the  computer  language  FORTRAN  IV.  The  analyzation  of  the  decision 
process  was  relatively  simple,  not  so  the  translation  into  FORTRAN  IV. 


Let  us  first  analyze  the  process  of  contour  suppressing  in  terms 
of  human  action,  the  interplay  between  the  memory  (brain  memory), 
the  optical  sensor  (eye),  and  the  mechanical  part  of  the  pen  movement 
(drawing  by  hand):  First  of  all  we  obtain  from  an  external  unit 
(somebody)  the  surface  information,  generally  in  terms  of  function  values 
at  grid  points,  together  with  the  boundary  outlines  and  the  information  in 
which  region  the  contours  should  be  drawn  (or  in  which  region  they  should 
be  deleted).  This  corresponds  to  the  information  represented  in  Fig.  2.8.4 
(+  surface  information). 


Figure  2.8.4  Boundary  outlines 

+  .  .  .  contour  region 
-  .  .  .  contour-free  region 
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1)  The  first  step  consists  of  a  global  "look"  at  the  situation 
in  order  to  get  an  impression  of  where  contours  have  to 
be  drawn.  If  no  other  means  is  available,  the  test-person 
will  assign  a  "YES"  to  all  contour  regions,  a  "NO"  or 
just  nothing  to  all  contour-free  regions,  and  a  "?"  to  the 
boundary  zone  which  is  neither  YES  nor  NO,  and  "store" 
this  YES-? -NO  pattern  in  his  brain.  Now  we  have  the 
following  information  available:  the  array  of  function 
values  at  the  grid  points  (  =  the  surface),  the  rough  pattern 
of  YES-? -NO  entries,  and  the  exact  boundary  outlines. 

2)  The  next  step  is  the  decision  where  we  should  start  search¬ 
ing  for  contours.  The  YES-? -NO  pattern  is  the  guideline 
for  this  decision.  Obviously,  one  does  not  start  amidst  a 
NO-  group  in  order  to  find  out,  after  the  calculation  of  the 
contour  line  element,  that  it  was  useless  anyway  because  of 
its  location  within  the  contour-free  region.  The  test-person 
would  probably  also  not  start  to  calculate  and  draw  contours 
amidst  a  YES-  group  --  this  would  speak  against  a  systematic 
solution  of  the  problem.  It  would  start  rowwise  {or  column¬ 
wise)  with  a  "?"  element,  an  element  from  a  boundary  zone. 
In  other  words,  the  test-person  starts  to  calculate  contours 

in  a  "zone"  along  the  boundaries. 

3)  He  calculates  the  line  element  of  the  contour  and  has  now 
to  intersect  this  line  element  with  the  boundary  (or  the 
boundaries).  For  this  purpose  he  needs  to  have  the 
"accurate"  information  about  the  boundary  (boundaries)  in 
th'is  small  region.  He  does  not  need  to  know  the  whole 
boundary  information,  but  only  a  "close-up"  of  a  very 
limited  part. 

4)  He  moves  his  eyes  closer  to  this  part,  concentrates  on  only 
two  lines,  the  contour  segment  and  the  small  part  of  the 


boundary,  and  forgets  for  some  time  all  other  regional  and 
global  information.  He  finds  the  intersection;  now  he  has 
to  decide  on  the  pen  position,  up  or  down.  Since  he  "forgot" 
in  the  meantime  all  global  information,  he  recalls  them 
again  by  moving  his  eyes  away  from  the  picture  in  order  to 
get  again  the  global  impression  of  the  YES-? -NO  pattern. 

This  pattern,  in  turn,  enables  him  to  choose  the  correct 
pen  position,  and  secondly,  to  make  him  trace  the  contour 
in  the  correct  direction  (the  direction  of  pen  down). 

5)  He  lowers  the  pen  at  the  point  of  intersection,  calculates  the 
next  contour  element;  if  this  happens  to  be  in  the  "?"  region 
he  repeats  steps  (3)  and  (4);  if  it  falls  in  the  YES  -  region 
he  continues  drawing  until  he  finds  the  next  11  ? 11  region. 

The  first  pen-up/pen-down  decision  is  sufficient  for  all  further 
ones  because  of  the  alternating  character  of  the  contour/ 
contour-free  region  sequence. 

These  5  steps  described  above  are  in  principle  "translated"  into 
FORTRAN  IV.  In  the  sequel  we  present  the  essence  of  this  translation. 

1)  As  described  in  Chapter  (2.7),  an  integer*  2  array  is 
assigned  to  the  array  of  bilinear  elements.  The  same 
array  is  used  for  storing  the  information  whether  a  particular 
element  (here  rectangle  on  which  a  bilinear  function  is  defined) 
happens  to  be  located  completely  inside  a  contour  region  (YES), 
completely  inside  a  contour-free  region  (NO),  or  if  the  element 
is  crossed  by  one  of  the  region  boundaries. 

The  boundary  information  consists  of  the  number  of 
boundaries,  the  information  whether  boundary  i  encloses 
a  contour  -  or  contour  free  region  (for  all  boundaries  i), 
the  boundary  coordinates  and  information  if  the  boundary 
coordinate  sequence  runs  clockwise  or  counter-clockwise. 
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In  the  sequel  all  elements  are  marked  by  a  negative 
number  (corresponding  to  the  question  mark  "  ? "  above), 
if  at  least  one  boundary  line  crosses  the  element.  The 
absolute  value  of  this  negative  number  serves  at  the  same 
time  as  counter  of  how  many  boundary  points  have  been 
found  (or  interpolated)  within  the  element  in  consideration. 
The  interpolation  of  boundary  points  is  necessary  if  their 
mutual  distance  is,  loosely  speaking,  too  large  relative 
to  the  size  of  the  element;  the  interpolation  is  linear,  for 
reasons  of  simplicity.  As  soon  as  the  first  boundary 
has  been  completed  (all  boundaries  need  to  be  closed)  and 
the  "border  zone"  has  been  marked,  the  region  is  filled 
up  with  "1"  's  (—  NO)  if  the  region  is  to  be  contour- 
free;  the  contour-regions  bear  a  "0".  After  all  regions 
have  been  processed,  a  pattern  of  "?",  "0",  "1"  is 

obtained.  This  pattern  represents  the  global  information 
of  where  contours  have  to  be  drawn  and  of  where  they 
should  be  deleted.  Fig.  2.8.5  shows  such  a  pattern 
associated  with  the  boundary  situation  of  Fig.  2.8.4. 

In  addition,  all  boundary  points  (primary  and  inter¬ 
polated)  which  will  be  used  to  calculate  contour-boundary 
intersections  are  stored  on  an  auxiliary  vector.  (Usually, 
not  all  boundary  points  are  needed,  unless  the  boundaries 
are  completely  inside  the  total  rectangular  plotting  range; 
e.g.  consider  the  case  of  the  complete  world  shore-outline 
data  bank  and  a  plot  restricted  to  Europe;  then  probably 
only  a  couple  of  thousand  out  of  the  total  80  000  boundary 
data  will  be  used. ) 

2)  The  next  step  is  already  a  part  of  the  actual  contouring. 

The  array  filled  up  with  negative  numbers  (boundary  zones), 
zeros  (contour  regions)  and  ones  (contour-free  regions)  is 
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Figure  2.8.5  Region  patterns 
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checked.  prior  to  the  function  value -contour  value  compari¬ 
son;  all  elements  filled  up  with  "1"  are  branched;  all 
contour  points  are  calculated  and  stored  for  successive 
"0"  elements. 

3)  As  soon  as  the  contour  touches  a  negative  valued  element 
(boundary  zone),  the  program  switches  to  a  "close-up" 
look  for  this  particular  area.  The  close-up  region  covers 
only  all  boundary  data  within  the  element  in  consideration 
plus  all  boundary  data  in  the  immediate  neighborhood  of 
this  element;  This  number  of  data  is  always  very  small 
indeed  (recall  that  one  element  has  a  maximum  size  of 
3x3  mm).  At  this  point  we  do  not  go  into  detail,  this 
would  be  a  little  bit  too  technical.  The  whole  close-up 
look  consists  in  a  series  of  logical  operations  involving 

a  tree  of  pointer  vectors;  the  only  calculations  are 
integer  additions  --  therefore,  this  procedure  is  extremely 
efficient  and  takes  very  little  computer  time. 

4)  Havii  issembled  all  local  boundary  points,  the  program 
continues  with  the  actual  calculation  of  intersections  which 
is  performed  in  a  stable  subroutine.  Such  a  line  inter¬ 
section  subroutine  was  available  to  the  author  in  form 

of  an  elegant  flowchart  (Neubauer,  1978).  The  subroutine 
following  this  flowchart  was  programmed.  If  no  intersection 
has  been  found,  the  program  continues  and  branches  to  the 
next  element.  If  one  intersection  has  been  found,  there 
is  the  problem  of  deciding  whether  at  this  particular  point 
the  pen  has  to  be  raised  or  lowered.  This  decision  is 
enabled  by  a  global  information  in  terms  of  a  logical  vector, 
assigned  to  the  boundaries,  which  provides  information  about 
the  region  behavior  in  a  direction  orthogonal  to  the  boundary 
(contour  or  contour-free);  e.g.,  in  Fig.  2.8.6  the  points  Pj 
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P2  are  contour  points,  #37,  ...,  #42  are  boundary  points 
of  the  boundary  B3;  they  represent,  for  this  particular 
element,  the  "close-up"  ;  (37,  42)  are  only  apart 

of  the  boundary  B3.  The  global  information  consists  in 
the  knowledge  of  which  side  of  the  boundary  B3  is  contour- 
free;  then  it  is  clear  that  at  S  the  pen  has  to  be  lowered 
and  the  connection  S-P2  has  to  be  plotted.  The  part  Pj-S 
of  the  contour  will  be  deleted. 


boundary  B3 

contour 


Fig.  2.8.6  Close-up  contour/simple  boundary 

The  case  of  a  single  intersection  point  per  element 
is  relatively  easy  to  handle.  Very  detailed  and  rough 
boundaries,  however,  may  cause  more  than  one  inter¬ 
section  points.  Moreover,  the  intersections  may  even 
refer  to  different  boundaries.  In  such  a  case  the  close-up, 
the  retrieval  of  the  logical  information  associated  with 
each  boundary,  the  calculation  of  all  intersection  points 
and  the  pen  position  logics  are  much  more  complicated. 


42  f' 
/ 
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boundary 

contour 

boundary 


B 
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Fig.  2.8.7  Close-up  contour/complicated  boundary 


E.g.:  In  Figure  2.8.7  the  points  Pj  and  F2  are 
contour  points  as  before.  But  now  there  are  two  boundaries 
involved,  B2  and  B3.  The  close-up  of  B2  is  identical  with 
B2  --  it  is  a  closed  boundary  within  a  single  element  con¬ 
sisting  of  the  points  #  76,  ...»  #80  (the  points  #76  and 
#  80  are  identical  because  the  boundary  is  closed).  The 
close-up  of  B3,  however,  shows  only  a  part  of  B3  --  the 
points  #116,  ...,  #125.  Altogether,  13  boundary  line 
elements  are  to  be  checked  for  an  intersection  with  the 
contour;  6  intersection  points  Sj  ...»  S&  have  been  found. 
Let  us  assume  that  boundary  B2  encloses  a  contour  region; 
clearly,  the  region  to  the  right  hand  side  of  boundary  B3 
has  to  be  also  a  contour  region  because  of  the  alternating 
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behavior  of  such  regions.  Let  us  now  intersect  all  boundary 

r 

line  elements  with  the  line  Pt  P2  ,  starting  with  boundary 
B2  in  increasing  order.  Then  it  can  be  seen  immediately 
that  the  sequence  of  intersection  points  (s3}  ,  i=l,  ...,  6 
is  more  or  less  arbitrarily  distributed  between  Pj  P2  (S5,  S4, 

S3>  Sj,  S2,  Sf,).  In  order  to  make  the  pen  position  decision 
it  is  necessary  to  know  just  one  intersection  point  together 
with  the  boundary  logics;  all  other  pen  positions  follow 
from  the  alternating  behavior.  This  alternating  continuation, 
however,  is  only  possible  if  the  actual  intersection  point 
sequence  is  known.  This  can  be  accomplished  by  first 
transforming  the  intersection  point  coordinates  into  a  new 
coordinate  system  whose  origin  is,  in  principle,  arbitrary 
and  whose  x-axis  is  parallel  to  the  line  Pj  P2,  and  then  re¬ 
ordering  the  intersection  points  according  to  increasing 
x-coordinates.  The  first  operation  is  just  a  rotation  of  the 
original  coordinate  system,  the  second  operation  is  not  so 
trivial  --it  involves  a  sorting  algorithm.  GSPP  uses  an 
algorithm  written  by  P.  Meissl  which  has  been  adapted  for 
this  particular  problem.  In  the  case  discussed  above  and 
graphically  shown  in  Fig.  2.8.7,  the  intersection  point  S5 
belonging  to  the  boundary  B3  would  have  the  smallest  trans¬ 
formed  x-coordinate  among  all  other  intersection  points. 

This  point  will  also  be  taken  for  the  pen-position  decision 
process:  From  the  global  information  it  is  known  that  the 
region  to  the  right  hand  side  of  the  boundary  B3  is  a  contour 
consequently,  the  segment  PtSs  has  to  be  plotted, 
deleted,  S3  plotted,  ...  and  so  on.  This  completes  step  (4), 

5)  In  th**  sequel  tthe  next  element  will  be  checked.  If  it  has  a 
value  m0”  assigned,  the  contouring  can  continue  without  the 
need  of  close-up's  and  intersection  calculations;  the  same 


is  true  for  values  "1"  --in  this  case  the  program  can  "jump" 
from  one  element  to  the  next  until  it  finds  all  elements 
with  a  negative  value  in  which  case  steps  (3)  and  (4)  have 
to  be  repeated. 

If  smooth  contours  are  requested,  the  program  has  to  "remember" 
all  intersection  points  including  at  least  the  pen  position  at  1  inter¬ 
section  point.  The  parts  of  the  contour  which  have  to  be  plotted  are 
then  smoothed  along  the  procedure  described  in  section  2.8.  1.  The 
following  Figure  2.8.8  serves  as  an  illustration.  Figure  2.8.9 
shows  that  even  for  regions  with  very  complicated  boundaries  the 
procedure  works  absolutely  correct  and  gives  total  resolution. 

2.8.6  Region  boundary  plot  .  Besides  the  suppressing  of  the 
contour  plot  within  certain  regions,  there  is  also  a  routine  built  in  which 
plots  the  actual  boundaries.  There  is  an  option  to  plot  some  of  them; 
there  is  the  option  to  assign  to  each  boundary  a  certain  linewidth. 

In  any  case  the  boundary  lines  will  be  clipped  off  at  the  border  of  the 
rectangular  plotting  domain.  The  boundary  points  will  be  connected 
linearly  throughout. 

The  clipping  process  has  been  designed  in  the  following  manner: 
assume  the  rectangular  plotting  domain  R  to  be  described  by  the  coordin¬ 
ates  [x0,Xj,  y0,  yt].  If  a  boundary  line  happens  to  be  completely  inside 
the  rectangle  R,  it  will  be  plotted  (after  an  optional  mapping)  immediately. 
If  a  boundary  line  happens  to  be  completely  outside  the  rectangle  R  (all 
x-coordinates  <  x0  or  all  x-coordinates  >  xt  and/or  all  y-coordinates  <  y0 
or  all  y-coordinates  >  yj),  the  program  skips  this  boundary  and  continues 
with  the  next  one.  If  none  of  the  above  conditions  is  true,  the  boundary 
either  crosses  the  border  of  the  rectangle  R  or  it  is  still  completely 
outside  R.  The  same  is  true  for  all  line  segments  of  the  boundary. 

The  following  number  pattern  is  assigned  to  the  nine  rectangular 
regions  (l  closed,  8  open)  generated  by  the  boundary  lines  of  the 
plotting  rectangle  and  its  infinite  continuations: 
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One  of  these  values  is  assigned  to  the  boundary  points  according  to 
its  position.  Then  it  can  be  immediately  verified  that  a  line  segment 
is  definitely  completely  outside  R  if  the  number,  which  results  from 
an  addition  of  the  codes  of  the  two  line  segment  points,  contains  at 
least  one  "2";  if  it  does  not,  it  may  be  completely  outside  but  might 
also  be  partly  inside  R,  and  therefore,  intersect  the  boundary  rectangle. 

If  the  sum  is  zero,  it  is  completely  inside  R.  E.g.: 

Pi  is  in  the  region  1010,  P2  in  0011;  the  code  sum 
is  1021  --  therefore,  the  line  Pj  P,  is  completely  outside 
R.  If  P2  is  in  0101,  then  the  code  sum  is  1111  -- 
therefore,  the  line  p^jcould  cross  the  boundary  of  R. 

If  the  line  segment  code  contains  a  "2",  the  line  segment  will 
be  deleted;  if  it  contains  only  zeros,  it  will  be  plotted;  if  it  contains 
at  least  one  "1"  but  no  "2",  the  program  finds  the  point(s)  of  inter¬ 
section  with  the  help  of  a  line  intersection  algorithm  described  in 
(Neubauer,  1978).  If  a  mapping  is  defined  and  requested,  all  boundary 
points  (and  intersection  points)  will  be  mapped;  their  connection  will  ‘ 

always  be  linear  regardless  of  the  mapping. 

j 

2.8.7  Plot  of  horizontal  and  vertical  axes.  The  axis  routine 
will  be  described  in  detail  in  Chapter  5  of  Fart  B.  Here  we  mention  that 
a  couple  of  axis  options  are  available  for  a  contour  plot. 
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Any  axis  plot  can  be  suppressed;  the  axis  can  be  obtained  on 
the  left  hand  side  and  on  the  bottom  of  the  contour  plot;  the  axis  can, 

* 

in  addition,  be  obtained  on  the  right  hand  side  and  on  the  top  of  the 
I  contour  plot.  The  axis  itself  consists  of  a  straight  line  (which  can 

also  be  suppressed),  tick  marks  and  scale  numbers.  The  tick  mark 
interval,  the  tick  mark  length  and  the  tick  mark  direction  can,  within 
certain  bounds,  be  chosen  by  the  user.  If  the  user  does  not  define 
a  tick  mark  interval,  it  is  taken  as  the  grid  interval  (for  the  surface 
representation).  The  scale  numbers  can  be  plotted  in  4  different 
directions  ( k #  ff/2  ,  k  =  0,  ...,  3);  the  symbol  heights  and  the 
number  of  decimal  places  can  be  defined  by  the  user,  otherwise  default 
values  will  be  assigned  (e.  g,,  such  that  at  least  3  significant  digits 
are  plotted).  A  scale  number  will  always  be  centered  with  respect 
to  its  corresponding  tick  mark,  where  the  actual  length  of  a  scale 
number  is  determined  by  equation  (2.8.2).  If  the  sequence  of  scale 
numbers  is  too  dense  (such  that  overlapping  would  occur),  certain 
scale  numbers  will  be  deleted.  In  addition  can  the  distance  between 
the  tick  marks  and  scale  numbers  be  chosen  by  the  user  (within 
certain  bounds).  Non-labeled  tick  marks  will  be  40%  shorter  than 
labeled  ones.  No  axis  will  be  plotted  if  a  contour  mapping  is  involved. 

2.8.8  Plot  of  a  grid  superimposed  on  the  contour  plot.  In 
order  to  make  graphical  interpolations  c  notion  values  between  tne 
plotted  contours  easier,  a  grid  plotting  routine  has  been  included  in 
GSPP.  The  following  grid  patterns  are  available  (Fig.  2,8.10): 
full  one,  dashed  line  with  arbitrary  dash  length  and  interval  between 

> 

dashes,  and  a  plot  of  open  crosses  at  the  grid  points.  The  grid 
t  distance  is  always  identical  with  the  tick  mark  interval.  If  a  grid  plot 

has  been  requested  and  no  grid  parameters  have  been  defined,  a  full 
line  grid  will  be  plotted.  No  grid  will  be  plotted  if  a  contour  mapping 
is  involved.  Further  details  about  a  grid  plot  can  be  found  in 
Chapter  6  of  Part  B. 
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2.8.9  Plot  of  data  superimposed  on  the  contour  plot. 

Especially  when,  a  surface  is  derived  from  a  set  of  homogeneous  data 
(noisy  or  not),  it  is  often  of  interest  to  plot  the  data  with  the  contours 
as  background;  this  enables  the  user  to  check  the  goodness  of  surface 
fit  relative  to  the  data.  But  also  a  plot  of  the  data  point  locations  and 
contours  might  be  quite  useful  if  one  wants  to  see  how  the  prediction 
works  in  data-free  areas  or  in  areas  with  poor  data  coverage.  In 
order  to  satisfy  all  needs,  GSPP  offers  the  following  3  options:  only 
data  position  symbol  plot,  data  position  symbol  plot  and  data  number 
plot,  data  position  symbol  plot  and  data  number  and  data  value  plot. 

The  position  symbol  can  be  chosen  among  all  symbols  available 
(differs  from  one  plotting  software  package  to  the  other)  --  some  of 
them  (0,  ....  13  for  IBM  software)  give  centered  symbols,  some  of 
them  dont.  The  data  number  and  the  data  value  will  always  be 
symmetrically  located  relative  to  the  data  position,  the  data  number 
above  the  data  position,  the  data  value  below  the  data  position.  The 
symbol  height  can,  within  a  certain  bound,  be  selected  by  the  user; 
the  same  is  true  for  the  number  of  decimal  digits  for  a  data  value 
plot.  In  any  case,  the  symbol  height  of  the  data  number  will  be  half 
the  symbol  height  of  the  data  value;  this-  is  to  distinguish  them  more 
easily. 

Besides  the  data  plot  there  is  also  the  possibility  to  plot  the 
predicted  data  in  the  same  way  as  described  above  (plot  the  grid  point 
locations  and  predicted  value  at  the  grid  points). 

If  a  mapping  is  involved,  the  data  coordinates  will  be  mapped 
in  the  same  way  as  the  contour  point  coordinates.  In  any  case  a  data 
point  outside  the  rectangular  plotting  domain  will  never  be  plotted. 
Figure  2.8.12  illustrates  a  data  plot  superimposed  on  a  contour  plot. 


I 


Figure  2.8.  12  Contour  and  data  plot 


2.8.10  Title  and  label  control  and  plot.  Titles  and  labels  are 


used  to  identify  plots.  For  this  reason  GSPP  offers  the  possibility  of 
consisting  of  up  to  10  title  lines,  as  well  as  single  line 
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horizontal  and  vertical  label  plots  (a  line  is  to  be  understood  as  a 
maximum  of  80  characters  =  information  stored  on  one  punched  card). 

The  height  and  line  width  of  the  symbols  can  be  chosen  by  the  user  within 
certain  bounds. 

The  title  plot  routine  itself  is  intelligent.  It  offers- the  title  to 
be  plotted  in  the  following  four  ways:  same  as  on  the  input  cards 
(default),  left  justified,  centered,  right  justified  regardless  of  how 
the  title  information  was  actually  punched.  If  the  title  length  exceeds 
the  actual  length  of  the  contour  plot,  the  symbol  size  will  be  reduced. 

If  the  reduced  symbol  height  happens  to  be  smaller  than  the  allowed 
lower  limit,  the  lower  limit  will  be  chosen  and  the  title  will  start 
on  the  left  side  regardless  of  the  requested  mode.  The  symbol  size 
will  also  be  reduced  automatically  when  the  total  title  height  exceeds 
upper  limits,  fixed  by  GSPP.  The  title  will  always  be  plotted  above 
the  contour  map;  its  position,  relative  to  the  upper  left  comer  point 
of  the  rectangular  plotting  domain,  can  be  chosen  by  the  user  --  if  not, 
default  values  will  be  assigned. 

Two  single  line  labels,  one  along  the  horizontal,  one  along  the 
vertical  axis  can  also  be  plotted.  As  long  as  the  label  length  does  not 
exceed  the  corresponding  length  of  the  contour  plot,  the  label  will  be 
centered  relative  to  the  contour  plot.  If  it  exceeds  this  length,  the 
symbol  height  will  be  reduced  automatically,  but  not  below  a  minimum 
height  defined  in  GSPP.  If  the  label  length,  with  minimum  symbol 
height,  is  still  bigger  than  the  corresponding  length  of  the  contour 
plot,  it  will  no  longer  be  centered;  instead,  it  will  be  left-justified. 

Usually,  the  labels  will  be  plotted  along  the  left  vertical  and  along 
the  lower  horizontal  axis;  however,  when  axes  are  to  be  plotted  all 
around,  also  labels  will  be  plotted  all  around.  In  this  case  the  title  will 
also  be  shifted  automatically  in  order  to  avoid  a  possible  overlapping. 

No  title  and  no  labels  will  be  plotted  if  a  contour  mapping  is  requested. 


LflT 
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2.8.11  Contours  of  surface  derivatives.  The  contouring  part 
of  GSPP  has  been  designed  primarily  for  a  contour  plot  of  the  surface 
itself.  The  surface  is  defined  as  a  bicubic  spline  function  which  is 
determined  by  function  values  at  the  grid  points  of  a  regular  rectangu¬ 
lar  grid.  For  the  contouring  it  is  of  no  concern  whether  the  grid 
values  were  known  in  advance  or  if  they  were  predicted  on  the  basis 
of  other,  probably  irregularly  distributed  inhomogeneous  data,  or  if 
they  are  function  values  derived  from  a  least- squares  regression 
polynomial.  All  of  them  represent  a  surface  which  is  interpreted  as 
a  bicubic  spline  function. 

Apart  from  surface  contours,  GSPP  also  offers  contours  of 
surface  derivatives  up  to  and  including  a  second  derivative  in  both 
coordinate  directions.  Therefore,  contour  plots  of  the  following  surface 
derivatives  can  be  requested  from  GSPP: 
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The  derivatives  are  calculated  according  to  the  set  of  equations  (2.6-7b), 
which  are  derivatives  of  a  bicubic  polynomial.  Since  a  bicubic  poly¬ 
nomial  is  twice  continuously  differentiable  with  respect  to  both  independent 

variables,  the  highest  derivative  offered,  D  f,  is  still  a  continuous 

xx  yy 

function;  it  is  a  continuous  and  piecewise  bilinear  function  and 

therefore,  a  hyperbolic  paraboloid  (cf.  Section  2.7).  All  lower  order 

derivatives  will  be  functions  of  higher  order,  and  therefore,  smoother. 

The  following  Figures  2.8.14(a,  b,  c)  show  a  contour  plot  of 

a  surface,  its  first  derivative  in  x-direction  and  its  highest  allowed 

derivative  D  f.  From  the  D  f  plot  one  can  see  very  clearly  the 
xx yy  xxyy 

contour's  tendency  to  run  indirections  along  the  asymptotic  lines  of 
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the  hyperbolas,  which  are  parallel  to  the  coordinate  lines,  and  along 
the  axes  of  the  hyperbolas  which  span  an  angle  of  45°  and  135°, 
respectively,  with  the  coordinate  lines.  These  features  are  particularly 
pronounced  in  Fig.  2.8.14c  because  of  the  big  size  of  the  grid  ele¬ 
ments  used  (2x2  cm).  This  happens  always  when  bilinear  elements 
are  contoured.  Since  the  bicubic  elements  are  approximated  by 
bilinear  elements,  the  same  is  true  even  in  the  contouring  of  surfaces. 
However,  the  size  of  the  subgrid  is  kept  so  small  that  the  effect  cannot 
be  seen  anymore. 


I 


l 
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3.  PROFILES 

The  term  "profile"  is  usually  thought  of  as  a  curve  which  results 
from  an  intersection  between  a  surface  and  a  vertical  plane.  Here  we 
understand  by  profile  a  curve  which  can,  but  needs  not  necessarily  be, 
a  curve  of  intersection  in  the  usual  sense.  If  a  curve  is  defined  by 
a  vector  of  pairs  (x,  y)  with  yj  =  y(x^)  function  values  at  x.  ,  we 
speak  about  an  "explicitly  defined"  profile.  If  a  prof  ile  is  defined  by 
a  surface  (or  data  which  are  to  represent  a  surface)  together  with 
profile  start  and  end  point,  we  speak  about  an  "implicitly  defined" 
profile.  If  more  than  one  explicitly  defined  profile  is  to  be  plotted 
in  the  same  frame,  we  will  speak  about  a  "multiple  profile". 

GSPP  can  handle  these  three  types  of  profiles  with  a  number 
of  options  and  additional  features  like  profiles  of  surface  derivatives, 
profile  derivatives,  profile  information,  and  many  others,  fully  auto¬ 
matically.  The  following  chapters  describe  these  features  in  detail. 
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3.  1  Explicitly  defined  profiles 

As  already  stated  in  the  foregoing  introduction,  an  explicitly 
defined  profile  is  thought  of  as  a  curve  which  is  sampled  by  a  vector 
of  increasing  arguments  (x.},  i  =  1,  I  together  with  a  vector 

of  function  values  (y\3i  i  =  1,  . . . ,  I.  Of  course,  there  is  a  number 
of  ways  to  connect  these  sampled  curve  points.  One  could  think  of 
a  polynomial  interpolation  or  least-squares  interpolation,  etc.  GSPP 
essentially  offers  three  kinds:  a  simple  sample  point  plot,  a  piecewise 
linear  interpolation,  and  a  cubic  spline  interpolation. 

Profile  derivatives  can  be  requested  up  to  the  second  order; 
derivatives  will  be  derived  from  a  cubic  spline  representation  of  the 
profile.  The  corresponding  routines  have  been  described  in  Sections 
2.6.1  and  2.6.3.  As  far  as  the  spectral  content  of  the  spline 
representation  is  concerned,  the  reader  may  consult  Section  2.5.1. 

The  actual  plot  of  the  profile  is  performed  within  a  window  in 
x-  and>  y-direction  (  argument  window  and  function  value  window).  If 
a  window  has  been  defined  by  the  user,  the  curve  will  be  clipped  off 

when  it  leaves  the  window.  The  clipping  algorithm  is  essentially  identical 
with  the  one  described  in  Section  2.8.6.  If  no  window  has  been 
defined,  the  minima  and  maxima  of  the  {x.}  and  (y.)  vectors  are 
assumed  to  coincide  with  the  bounds  of  the  window. 

3.1.1  Optional  profile  procedures.  The  profile  interpolation 
can  be  done  piecewise  linearly  or  by  an  interpolating  cubic  spline. 

For  reasons  of  stability,  the  spline  is  calculated  piecewise  and  over¬ 
lapping  if  the  number  of  data  points  is  too  large  (see  Section  2.8.1). 

If  no  interpolation  is  requested,  the  profile  points  will  be 
marked  by  a  centered  symbol  which  the  user  can  choose. 

The  plot  of  horizontal  and  vertical  axes  is  identical  to  that  one 
described  in  Section  2.8.7  with  the  exception  that  no  horizontal  axis 
can  be  plotted  above  the  profile  plot. 
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Analogous  to  the  superimposed  grid  in  the  contour  plot,  small 
crosses  can  be  plotted  at  the  intersections  of  grid  lines  (horizontal 
and  vertical  lines  at  tick  mark  intervals).  These  crosses  will  only 
be  plotted  below  the  profile.  In  addition,  there  is  the  option  to 
draw  a  horizontal  "zero-line"  whenever  the  plotting  window  for  function 
values  contains  the  zero  point. 

Title  and  label  plots  are  also  identical  to  those  described 
in  Section  2. 8.  10, 

As  far  as  the  plot  of  profile  derivatives  is  concerned,  there 
is  the  option  to  plot  the  first  or  the  second  derivative.  Since  the 
derivatives  are  taken  from  a  cubic  spline  representation,  the  second 
derivative  will  still  be  continuous  --  it  is  a  piecewise  linear  function. 
The  following  Figures  3.  1.  1  (a,b,  c)  serve  as  illustration  examples. 
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3.2  Multiple  profiles 

In  the  case  that  one  or  more  profiles  have  to  be  compared  with 

each  other,  it  is  advantageous  to  plot  all  of  them  in  the  same  frame. 

In  that  case  there  is  one  natural  restriction:  the  plotting  window  for 

all  profiles  have  to  be  identical  and  to  be  known  before  the  first 

profile  is  plotted.  If  no  such  plotting  window  has  been  defined,  GSPP 

assumes  the  minima  and  maxima  of  the  first  {x.}  and  Cy.]  vectors  to 

1  i 

represent  the  window.  However,  in  this  case  an  unwanted  clipping  of 
further  profiles  (in  the  same  frame)  occurs  whenever  it  exceeds  the 
x-  and  y-  minima/maxima  of  the  first  profile. 

In  order  to  simplify  the  identification  of  the  different  profiles, 
two  additional  features  (relative  to  the  single  profile  plot)  have  been 
built  in:  the  profiles  can  be  marked  by  centered  symbols,  each  pro¬ 
file  by  a  different  one;  the  number  of  symbols  relative  to  the  number 
of  profile  data  points  can  be  selected  by  the  user.  Furthermore,  there 
will  be  a  legend  plotted  at  the  right  end  of  the  profile  which  consists 
of  a  list  of  symbol-input  sequence  identifications.  The  symbol  plot 
along  the  profiles  and  the  corresponding  identification  plot  can  be 
suppressed  if  so  desired.  Figures  3.2.1a,  b  show  multiple  profile 
plots  without  and  with  symbol  suppressing. 

The  linewidth  can  be  changed  from  profile  to  profile.  Also 
the  derivative  can  be  changed  from  profile  to  profile. 

All  other  optioned,  procedures  like  axes  plots,  title  and  label 
plots,  etc.,  are  identical  to  the  single  profile  case. 

3. 3  Implicitly  defined  profiles 

A  profile  derived  from  a  surface  is  herein  called  an  implicitly 
defined  profile.  The  surface  is  assumed  to  be  a  bicubic  spline  surface 
defined  on  a  regular  rectangular  grid.  If  the  data  are  not  regularly 
distributed,  the  spline  surface  will  be  predicted  first.  Of  course,  this 
would  not  be  necessary  if  only  a  single  profile  is  calculated;  usually 
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one  needs  a  number  of  them  and  this  is  why  we  have  decided  to 
predict  the  whole  surface.  The  prediction  processes  are  identical 
with  the  ones  described  in  Section  2.2,  the  spline  representation  is 
described  in  Section  2.4.2.  The  other  essential  information  is  the 
start-  and  end-  point  coordinates  of  the  profile  in  consideration. 

The  profile  needs  not  necessarily  be  a  surface  profile  -- 
it  can  be  a  profile  of  any  derivative  of  the  surface 


3gf(x,  y) 
3*«»ay«*  ’ 

or  also  a  profile  of 


daf(x,  y) 

*  a 


8 S  4  *  »  Otz  s  2  ,  a  -  +  ctz 


(3.3-1) 


a  *  2 


(3.3-2) 


with  ds  the  line  element  of  the  straight  line  connecting  start-  and  endpoint 

of  the  profile.  AL1  derivatives  (3.3-1)  are  discussed  in  Section  2. 6. 3; 

the  derivatives  (3.3-2)  are  simply  the  projection  of  the  surface  gradient 

onto  the  unit  vector  e  with  direction  Pt-P2  (Pt. . .  start  point,  P2... 

T 

endpoint),  e  =  (cos  A,  sin  A) 


of.  _ 
d  s 

and  for  the 


7f*e  =  f  cos  A  +  f  sin 


x  y 

second  order  derivative  we 


T 

e 


Fe 


A, 

obcain 


(3.3-3) 


(3.3-4) 


with  the  second  order  tensor 


then  (3.3-4)  has  the  form 
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32f 


f  cos2  A  +  2f  cos  A  sin  A  +  f  sin2 A  . 
xx  xy  yy 


(3.3-4)  ' 


Since  the  bicubic  spline  is  twice  continuously  differentiable  with  respect 
to  x  and  y,  the  derivative  3f/6s  is  continuous  and  has  even  a  continuous 
derivative;  32f/^s2  is  only  continuous. 


The  actual  calculation  of  the  profile  is  done  pointwise  along 
the  line  Pj-P2  with  a  sampling  rate  of  2  points/mm;  the  profile  is 
then  a  linear  connection  of  all  profile  points. 

3.3.1  Implicit  profile  procedures.  There  are  a  couple  of 
features  which  are  only  connected  with  the  implicit  profile  plots.  First 
of  all,  the  request  for  a  smooth  profile  point  interpolation  will  be 
ignored  because  the  function  itself  is  already  smooth  and  the  profile 
point  sampling  rate  is  sufficiently  high;  therefore,  a  smoothing  would 
show  up  no  difference  to  the  linear  point  connection. 

Since  GSPP  is  capable  of  calculating  and  plotting  100  different 
profiles  by  just  a  single  call,  it  is  absolutely  necessary  to  identify 
the  different  profiles.  This  is  done  automatically  by  a  start-  and 
endpoint  message  which  appears  at  the  bottom  of  the  profile  (see  Figs. 
3.3,  1  a,  b).  This  message  will  always  appear  and  will  be 

centered  unless  its  length  exceeds  the  profile  length  even  with  minimal 
symbol  height  --in  this  case  it  will  be  plotted  left  justified. 

If  the  coordinates  of  the  start-  and/or  endpoint  are  such  that 
the  profile  happens  to  be  outside  the  actual  surface  domain  by  107<> 

(this  is  the  rectangular  region  on  which  the  bicubic  spline  surface  is 
defined),  the  plot  will  not  be  executed  for  this  particular  profile  and 
an  informative  message  will  be  edited  on  the  line  printer  (or  any  other 
selected  output  unit). 

Another  difference  fcr  the  explicit  profile  plots  is  the  scaling 
of  the  horizontal  axis.  There  are  essentially  two  kinds  of  scaling 
offered  by  GSPP:  a)  the  scaling  is  done  according  to  the  actual 
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distance  Pj-P2  with  a  tick  mark  interval  which  can  be  either  chosen 
by  the  user  or  can  be  determined  automatically  ; 

b)  a  value  1.0  is  assigned  to  the  total  length  Pj-P2  and  the  scaling  is 
done  accordingly. 

A  last  difference  is  the  grid  plot:  recall  that  a  rectangular 
plotting  window  can  be  selected  in  the  contour  plot  (see  Section  2.8.4). 
When  such  a  window  has  been  defined  and  the  profile  runs  within  the 
window,  crosses  will  be  plotted  at  the  intersections  of  grid  lines 
(horizontal  and  vertical  lines  at  tick  mark  interval).  At  the  point 
where  the  profile  crosses  a  window,  a  vertical  bar  will  be  plotted 
with  the  same  linewidth  as  the  profile  itself;  outside  the  window  but 
inside  the  surface  domain  tiny  crosses  will  be  plotted  instead  of 
normal  size  crosses  in  order  to  indicate  that  this  part  of  the  profile 
runs  within  an  area  for  which  no  contour  plot  has  been  performed. 

No  profile  and  no  crosses  will  be  plotted  if  the  profile  leaves  the 
surface  domain  up  to  the  above  mentioned  10%  limit.  Crosses  will 
only  be  plotted  below  the  profile.  The  following  Figures  3.  31(  a,  b,  c) 
illustrate  different  kinds  of  profiles. 

4.  THREE-DIMENSIONAL  SURFACE  REPRESENTATIONS 

A  3-D  representation  of  a  surface  is,  in  this  context,  to  be 
understood  as  a  projection  of  a  two-dimensional  surface  which  is 
embedded  in  a  three-dimensional  Euclidean  space,  onto  a  plane.  The 
plane  can  be  arbitrarily  oriented  in  space.  The  surface  can  be  given 
either  explicitly  by  function  values  at  the  grid  points  of  a  regular 
rectangular  grid  or  by  irregularly  distributed  data  in  which  case  a 
prediction  algorithm  (Section  2.2)  takes  care  of  the  prediction  of 
function  values  at  all  the  grid  points.  Again,  the  surface  is  considered 
to  be  a  bicubic  spline  surface  defined  by  the  function  values  at  the 
grid  points  (see  Section  2.6.2).  The  bicubic  spline  surface  is  then 
approximated  by  small  bilinear  elements;  the  function  values  at  the 
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Figure  3,  3.  lc:  Profile  of  the  surface  derivative  in  profile  direction. 


subgrid  points  (the  grid  on  which  the  bilinear  elements  are  defined) 
are  interpolated  spline- surface  values.  Up  to  this  step  the  3-D  pro¬ 
cedure  does  not  differ  from  the  contouring  described  in  Chapter  2. 

In  the  contouring  procedure,  these  bilinear  elements  have  to  be 
intersected  with  a  number  of  horizontal  planes;  in  3-D  representations 
the  bilinear  elements  are  projected  onto  a  plane  with  arbitrary- 
orientation  in  3-D  space. 

4.  1  The  projection 

In  principle,  one  could  think  of  any  kind  of  projection  of  the 
surface  onto  the  plane.  GSPP  assumes  that  the  projection  is  an 
axonometric  projection,  where  the  center  of  projection  is  located  at 
infinity.  Consequently,  two  parallels  will  remain  parallel  after  the 
projection.  (The  projection  equations  used  in  GSPP  could  easily 


be  changed;  therefore,  it  is,  in  principle, possible  to  obtain  two 
perspective  projections  from  two  different  centers  of  projection 
plotted  in  two  complementary  colors  [red  and  green],  which,  viewed 
with  anaglyph  glasses,  give  a  three-dimensional  impression  of  the 
surface. ) 

The  axonometric  projection  used  in  GSPP  is  defined  in  the 
following  way:  Assume  the  surface  z  =  z(x,  y)  to  be  defined  on  a 
rectangular  domain  D,  with  the  lower  left  point  of  the  rectangle 
coinciding  with  the  origin  (x=0,  y=0),  and  the  sides  of  the  rectangle 
parallel  to  the  coordinate  lines.  Let  a  plane  pass  through  this  origin. 
The  orientation  of  the  plane  is  defined  by  two  angles,  the  longitude  X 
and  the  co-latitude  9  in  the  following  way:  a  coordinate  system 
(x,  y,  z)  is  associated  with  the  plane  with  x  =x,  y  =  y,  z  =  z  if 
X  =  0  and  9=0.  Let  now  the  surface  be  fixed  and  let  the  plane  rotate 
around  the  z-axis  by  the  angle  X  in  the  positive  direction  (looked  upon 
from  the  origin,  this  is  a  clockwise  rotation;  looked  upon  from  a 
point  above  the  surface,  this  is  a  counter-clockwise  rotation).  The 
rotated  coordinate  system  will  be  called  (x1,  y',  z')- system.  Any  point 
Pwith  coordinates  (x,  y,  z)  will  have  coordinates 


The  second  rotation  will  be  performed  around  the  x'-axis  by  an  angle  9 
in  the  positive  direction  (as  defined  above).  The  so  obtained  coordinate 
system  is  called  (x,  y,  z)-  system.  In  order  to  be  clear,  after  the 
transformation  , the  two  planes  z=0  and  z  =  0  span  an  angle  9  with  each 
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other.  Any  point  P  with  coordinates  (x,  y,  z)  will,  therefore,  have 
coordinates 


and  written  explicitly  , 

x  =  x  cos  X  +  y  sin  X 

y  =  -xcos0  sin  X  +  ycos  0  cos  X  +  z  sin  0  (4.1-1) 

z  -  x  sin  0  sin  X  -  y  sin  9  cos  X  +  z  cos  0  . 

The  next  step  is  the  projection  of  the  point  P,  from  the  center  of 
projection  on  the  z-axis  at  infinity,  onto  the  (x,  y)  -plane,  which  gives 
the  Cartesian  coordinates  of  the  image  point  P*  of  P, 


(4.1-2) 


Combining  (4,1-1)  and  (4.1-2),  replacing  the  co-latitude  0  by  the  latitude 
9  =  90  -  0  ,  allowing  a  coordinate  shift  in  the  (x*,  y*)  -  system  and  scale 
factors  cj  for  (x,  y)  and  c2  for  z  we  obtain  the  final  projection  equations 
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x*  =  x0  +  (x  cos  X  +  y  sin  X)  ct 

#  * 

y  =  y0  +  (-x  sin  9  sin  X  +  y  sin  9  cos  X)  ct  +  z  cos  9  c2  . 

Another  interpretation  of  (4.  1-3),  apart  from  the  shift  (xj,  y£  )  and 
the  scale  factors,  is  the  following:  The  surface  domain  is  a  rectangle 
in  the  equatorial  plane  with  the  x-axis  through  the  origin  and  with 
orientation  angles  9  =  0  and  X  =  90°  ;  the  y-axis  passes  also  through 
the  origin  with  orientation  angles  9  =  0  and  X  =  180°;  then  the  sur¬ 
face  is,  according  to  equation  (4.  1-3),  looked  upon  from  a  point  at 
infinity  with  coordinates  9  and  X  on  the  unit  sphere. 

4. 2  Scale  and  shift 

In  a  3-D  plot  the  user  usually  faces  the  difficulty  of  reducing 
the  coordinates  x,  y  and  the  corresponding  function  values  z  such  that 
the  figure  "looks  nice".  Moreover,  there  is  the  problem  of  getting  the 
plot  on  a  particular  place  of  the  plotting  sheet.  This  would  make  a 
number  of  decisions  and  calculations  necessary  if  the  data  are  already 
regularly  distributed  on  the  rectangular  grid  and  if  they  are  known  to 
the  user.  He  would  first  have  to  find  the  minimum  and  maximum  of 
the  function  values,  then  make  the  necessary  projections  in  order  to 
find  out  the  plotted  size  of  the  surface  and  so  forth.  If  the  data  are 
irregularly  distributed,  probably  heterogeneous  and  noisy,  such  a 
decision  usually  becomes  a  pure  guess.  GSPP  is  smart  enough  to  do 
this  job  if  the  user  wants  it  to  be  done,  moreover,  he  can  still  make 
all  or  part  of  his  decisions  --  GSPP  will  accept  them  if  they  are 
consistent  and  reasonble,  and  will  reject  them  and  replace  them  by 
reasonable  ones  if  they  were  unreasonable. 

Let  us  briefly  describe  how  GSPP  finds  scales  and  shift  parameters 
The  scale  factor  cj  is  determined  such  that  a  2-D  plot  would  have  an 
across  the  plot  sheet  extension  of  10  cm.  Since  an  axonometric  pro¬ 
jection  reduces  lengths  (or  keeps  it  constant),  a  square  of  10  x  10  cm 
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never  exceeds  an  along  or  acioss  maximum  of  20  cm.  The  next 
step  consists  in  the  calculation  of  the  minima  and  maxima  of  the  surface 
function  values.  These  extreme  values  are  projected  using  a  simpli¬ 
fied  form  of  (4. 1-3)  and  the  scale  factor  c2  for  the  function  values  is 
determined  such  that  the  projected  difference  of  these  extrema  is 
7.5  cm.  If  the  values  (0,0)  are  assigned  to  the  direction  of  view 
(X,9  ),  GSPP  will  interpret  them  as  "no  input"  values  and  will  assign 
default  values  X  =  -60°,  9  =  30°  which  usually  gives  a  nice  view  of  the 
surface,  hi  order  to  find  the  actual  along  and  across  (the  plot  sheet) 
extensions,  it  is  necessary  to  find  the  projected  coordinates  of  the 
four  corner  points  of  the  rectangle  with  minimal  and  maximal  function 
values  assigned  to  them.  This  gives  obviously  only  upper  bounds  of 
the  3-D  plot  size  which  are  then  used  to  find  the  appropriate  shift 
parameters  (xj  ,  y£  ).  (The  determination  of  the  exact  extension 
across  the  plot  sheet  would  make  the  projection  of  all  surface  points 
necessary.  This  is  a  rather  time-consuming  task  and  should,  therefore, 
be  avoided. ) 

4.3  The  3-D  plot 

As  soon  as  the  bicubic  spline  surface  is  available,  a  piecewise 
bilinear  approximation  will  be  calculated  by  a  simple  interpolation  pro¬ 
cedure.  The  size  of  the  bilinear  elements  can  either  be  chosen  by  the 
user  or,  otherwise,  will  be  selected  by  GSPP.  The  interpolation  is 
done  first  rowwise,  then  columnwise:  the  first  along  profile,  a  con¬ 
tinuous  and  piecewise  linear  function  is  interpolated  from  the  surface, 
is  projected  using  (4.  1-3)  and  plotted;  then  the  next  parallel  profile  is 
interpolated  and  projected.  This  profile  and  further  profiles  have  to 
pass  a  hidden  line  algorithm  which  determines  the  actual  visibility  of 
the  profile;  those  profiles  or  portions  of  profiles  which  are  hidden  by 
previous  profiles  are  masked  and  will  not  be  plotted.  The  hidden  line 
algorithm  used  in  GSPP  is  a  modified  version  of  that  one  described  in 
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(Watkins,  1973).  In  this  way  all  horizontal  and  vertical  profiles  of 
the  surface  are  treated.  After  completion  a  frame  is  plotted  which 
also  passes  the  hidden  line  algorithm.  The  frame  plot  can  be  suppressed. 

4. 4  Additional  and  optional  procedures 

In  any  case  ,  information  about  the  viewing  direction  will  be 
plotted  in  the  lower  left  comer  of  the  plot.  This  information  gives  the 
used  values  of  the  longitude  and  latitude  of  the  direction  of  view  in 
degrees  and  minutes. 

A  rectangle  will  be  drawn  around  the  entire  plot;  it  can 
be  suppressed. 

At  the  top  a  title  consisting  of  no  more  than  10  title  lines  can 
be  plotted  in  the  same  way  as  for  the  contour  and  profile  plots  (with 
options;  left  justified,  centered,  right  justified,  or  as  on  the  input  cards). 
The  following  Figures  4. 4.  1  (a,  b)  show  a  contour  plot  and  a  corresponding 
3-dimensional  view. 


Figure  4.  4.  la: 
Contour  plot 
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It  could  have  been  anticipated  at  the  beginning  of  this  report 
(and  actually  shown  in  PART  A)  that  contouring,  profile  finding  and 
drawing,  and  the  generation  of  a  3-D  plot  consists  of  many  pro¬ 
cedures#  The  user  of  GSPP  need  not  ba  confronted  with  the  solution 
of  these  problems;  however,  many  of  these  procedures  (or  sub¬ 
routines)  cam  be  used  in  an  isolated  form  and  might  be  quite  useful 
as  such  or  as  parts  of  other  algorithms.  Therefore,  this  section  of 
the  report  discusses  all  the  procedures  and  algorithms  which  are 
independent  but  integrated  parts  of  GSPP. 

1.  Organization  of  data 

As  part  of  the  contouring  algorithm,  the  data  sorting  and 
retrieving  procedure  was  briefly  described  in  Section  2.1  of  part  A.  By  data 
organization  we  under  stand  an  algorithm  which  generates  pointer 
vectors  based  on  the  two-dimensional  distribution  of  the  data.  These 
pointer  vectors  should  make  a  very  fast  retrieving  of  data  within  a 
specified  array  possible.  The  generation  of  the  pointer  vectors  should 
take  as  little  time  as  possible  since  a  huge  number  ov  data  will  be 
organized.  The  following  is  a  description  of  the  method  used  by  GSPP. 
The  organization  is  performed  by  the  subroutine  OAF. 

First,  the  working  (or  organization)  domain,  which  is  assumed 
to  be  a  rectangle,  has  to  be  defined  by  its  lower  and  upper  x-  and  y- 
coordinates  (the  sides  of  the  rectangle  are  assumed  to  be  parallel  to 
the  coordinate  lines).  This  rectangle  is  divided  into  M*N  subrectangles 
of  equal  area  (M  in  x-direction,  N  in  y-direction).  The  program  finds 
for  each  coordinate  pair  (x,  y)  the  corresponding  element  (m,  n).  After 
a  couple  of  operations  (mainly  integer  additions  and  subtractions)  four 
pointer  or  counter  vectors  are  generated,  IC1(.),  ...,  IC4(,): 

The  vector  IC1(.)  has  a  length  equal  to  the  number  of  data 
and  contains,  after  completion,  the  element  index  i  corresponding  to 
each  data.  The  element  index  i  results  from  the  two  subrectangle 
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indices  (m,  n),  associated  with  a  data  point  (x,  y),  and  is  calculated 
according  to 

i  =  (m-  1)N  +  n, 

which  means  i  increases  rowwise  if  x  is  oriented  northwards,  i  increases 
columnwise  if  x  is  oriented  eastwards.  Therefore,  i  =  ICl(j)  is  the 

index  of  the  subrectangle  in  which  the  data  point  (x.,  y.)  is  located. 

3  3 

If  the  data  happens  to  be  outside  the  working  rectangle,  an  index 
M*N  +  1  will  be  assigned. 

The  vector  IC3{.)  is  a  counter  vector  of  length  M*N+1;  IC3(i) 
equals  the  number  of  data  in  the  subrectangle  with  index  i.  IC3(M#N+1) 
is  equal  to  the  number  of  data  outside  the  rectangular  working  domain. 

The  vector  IC4(. )  is  an  auxiliary  vector  of  length  M#N+1  ;  its 
elements  are  partial  sums  of  the  elements  of  IC3{.):  IC4(1)  =  1, 

IC4(k)  =  IC4(k- 1)  +  IC3(k- 1),  k  =  2,  ...,  M#N+1. 

Finally,  the  most  important  vector  is  IC2(. );  it  is  organized  in  such 
away  that  the  first  data  in  the  subrectangle  with  index  i  has  the  original 
index  IC2(lC4(i));  its  length  is  equal  to  the  number  of  data. 

The  data  retrieving  process  runs  then  as  follows.  Assume  one 
wants  to  know  all  data  which  are  located  within  the  element  i:  there 
are  altogether  IC3(i)  data  in  this  element;  the  index  of  the  first  data 
is  IC2(lC4(i)),  the  index  of  the  second  data  is  IC2(lC4(i)+l)  ,  and  so 
forth;  the  last  data  has  the  index  IC2(lC4(i)  +  IC3(i)-l).  If  there  are 
no  data  within  the  element  i,  then  IC3(i)  =  0  and  the  index  counting 
would  be  one  step  backward.  Therefore,  whenever  there  are  no  data 
within  a  particular  element,  the  index  retrieval  described  above  does 
not  apply. 

The  following  example  may  illustrate  the  foregoing  (Table  1.  1 
and  Fig.  1. 1).  This  kind  of  data  organization  is  extremely  fast  because 
there  are  only  very  simple  and  very  few  operations  involved;  the  data 
itself  are  not  shifted  --  they  remain  on  their  original  storage  locations. 


non  on  n—  ooooooo  ooo 
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A  DATA  ORGAN IZAT I  ON  EXAMPLE 
SUBROUTINES  USED  :  OAF 
OUTPUT  UNIT  :  6 

COMMON  /XDAT/XC  30)  /YDAT/Y( 30)  /USORPR/XMI , XMA, YMI , YMA, DUMMY!  2) , M, 

*  N  /S0RT1/ICK30)  /SORT2/IC2(30)  /S0RT3/TC3(  13) 

*  /S0RT4/ IC4( 13) 

TEE  VECTORS  IC1( . ) ,  . . . ,  IC4( .)  NEED  NOT  HAVE  THEIR  PROPER  DIM. 

AS  LONG  AS  THE  NUMBER  OF  DATA  IS  LESS  THAN  OR  EQUAL  TO  1000  AND  AS 
LONG  AS  THE  NUMBER  OF  SORT  ELEMENTS  IS  LESS  TRAN  OR  EQUAL  TO  1000 

THE  COORDINATES  OF  THE  DATA  POINTS  HAVE  TO  BE  STORED  ON  XC . )  AND 
Y( . ) .  HERE  WE  GENERATE  THESE  DATA  : 

THE  NUMBER  OF  DATA  IS  NDAT 
NDAT=30 
DO  1  1=1,  NDAT 
X(  I)  =SIN(  I*2.)+0.1 
Y(  I)=C0S( IS3.)-0. 1 
DEFINE  THE  WORKING  RECTANGLE 
XMI=-1. 

XMA3 1 . 

YMI=-K 
YKA=  1 . 

DEFINE  THE  NUMBER  OF  SUBRECTANGLES  (SORT  ELEMENTS)  FOR  THE  DATA 
ORGANIZATION 
M=3 
N=4 

THE  NUMBER  OF  ELEMENTS  IS  12 

CALL  THE  DATA  ORGANIZATION  SUBROUTINE  OAF  (ONE  ARGUMENT  =  NUMBER 
OF  DATA) 

CALL  OAF(NDAT) 

C  PRINT  THE  RESULTS  (INDEX,  DATA,  VECTORS  ICH.),  ....  IC4(.)  ) 
WRITE(6,60G0) 

NN1=M*N+1 
DO  2  1=1, NDAT 

IF< I.LE.NN1)  WRITE! 6, 6001)  I,X(  I) ,Y(  I) , IC1(  I) , IC2( I) , IC3( I) , IC4( I) 
IF(I.GT.MNl)  WRITE(6,6002)  I,X( I) , Y( I) , IC1( I) , IC2( I) 

2  CONTINUE 

6000  FORMAT(  1H0.4X,  *  I  XU)  Y(I)  ICK  I)  IC2<  I)  IC3(  I)  IC 

14(  I)  ’  ,//) 

6001  FORMAT! 1H  ,4X, I2,2(2X,F6,2) ,3X,4I8) 

6002  FORMATUH  ,4X.  I2,2(2X,F6.2)  ,3X,2I8) 

STOP 

END 


Program  for  Table  1.  1 


Figure  1.  1:  ’ 

Data  distribution 
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Table  1.  1  : 


Data  organization 
for  data  shown  in 
Fig.  1.  1  . 
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2.  Prediction 

The  basic  principles  of  the  prediction  methods  offered  by  GSPP 
have  been  shortly  described  in  section  2.2  of  PART  A.  In  this  chapter 
the  practical  aspects  of  least- squares  prediction  will  be  discussed. 

The  idea  of  least- squares  collocation  is  to  take  into  account  all 
gravity  field  information,  represented  in  terms  of  a  heterogeneous  data 
set,  for  the  prediction  of  other  gravity  field  quantities;  in  the  ultimate 
case  of  collocation,  model  parameters  can  be  incorporated  into  the 
solution. 

This  fine  concept  can  hardly  ever  be  fully  realized  in  practice; 
reality  demands  sacrifices.  The  following  facts  make  the  unified  solution 
a  prohibitive  task:  covariance  matrices  are,  in  contrast  to  network 
normal  equation  matrices,  full  matrices,  (in  collocation  one  has  to  deal 
with  the  continuum  "gravity  field",  in  network  problems  with  the  discrete¬ 
ness  of  a  continuum.)  The  size  of  the  matrix  depends,  again  in  contrast 
to  network  problems,  on  the  number  of  data.  The  more  data  we  have, 
the  better  we  have  sampled  the  gravity  field  and  -curiously  enough-  the 
more  difficult  it  becomes  to  determine  the  gravity  field:  the  instability 
of  the  covariance  matrix  increases  with  the  data  density  since  the  equations 
become  nearly  linear  dependent.  Apart  from  the  instability  there  is  the 
problem  to  store  the  matrix  -  a  very  serious  problem  if  more  than  a 
couple  of  thousand  data  are  involved,  not  to  speak  about  the  actual  in¬ 
version  or  calculation  of  the  solution  vector.  Last,  but  by  no  means 
least,  there  remains  the  actual  calculation  of  the  signal  (the  gravity 
field  quantity)  together  with  its  estimated  error  at  a  huge  number  of  grid 
points  --  we  have  to  keep  in  mind  what  we  actually  want:  the  determina¬ 
tion  of  a  gravity  field  surface  which  is  sufficiently  well  represented  by  an 
array  of  function  values  together  with  an  appropriate  interpolation  function. 

How  can  we.  overcome  these  problems,  what  are  the  consequences? 
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A  widely  applied  and  generally  accepted  practice  is  data  selection. 

(if  in  trouble,  select  --  in  analogy  to  Jeffreys  recommendation  "If  in 
doubt,  smooth".)  It  is  an  obvious  and  well  understood  fact  that  the  data 
in  the  immediate  neighborhood  of  the  prediction  point,  in  general,  contributes 
the  most  and  remote  data  very  little  to  the  prediction  of  the 

gravity  field  quantity.  (This  does  not  hold,  e.  g.  ,  for  the  prediction  of 
geoidal  heights  from  gravity  anomalies. )  Local  samples  are  considered 
rather  than  regional  or  global  ones.  The  number  of  data  used  for  a  pre¬ 
diction  depends  primarily  on  the  problem  (and  sometimes  on  the  personal 
taste).  Rapp  (1979)  uses  only  a  very  few  (around  5)  altimeter  data  for 
the  prediction  (or  rather  interpolation)  of  an  array  of  geoidal  heights  and 
in  (Rapp,  1978)  some  200  points  for  the  recovery  of  mean  gravity  anomal¬ 
ies  from  altimeter  data.  Schwarz  (1976)  stresses  the  fact  that  neighbor¬ 
hood  -  data  are  the  essential  information  and  suggests  a  data  selection. 
Lachapelle  (1977)  considers  some  one-to-two  hundred  gravity  anomalies 
and  deflections  of  the  vertical  for  the  combined  solution  collocation  and 
integral  formulas. 

The  consequences  are  as  follows:  local  collocation  solutions 
prevent  an  estimation  of  regional  and  global  parameters  like  datum  shifts 
for  obvious  reasons.  Collocation  in  the  local  mode  can  only  provide  a  less 
than  optimal  gravity  field  solution  since  an  optimal  soiutiai  would  require 
all  data  to  betaken  into  consideration.  This  is  the  price  we  pay  for  a 
gain  in  matrix  stability,  limitation  of  storage  requirements  and  for  keep¬ 
ing  the  computation  time  at  an  acceptable  level. 

In  gravity  field  surface  prediction,  the  factor  time  plays,  apart 
from  the  others  discussed  above,  a  particularly  important  role.  This 
is  why  only  local  solutions  can  be  envisioned  under  the  present  circum¬ 
stances. 

The  surface  prediction  algorithum  of  GSPP  is  designed  for  the 
local  mode  only.  It  considers  up  to  100  data  in  the  neighborhood  of  the 
prediction  point,  (This  number  has  been  kept  low  because  of  the 
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storage  limitations  of  tne  used  computer;  an  increase  to  200  or  even 
300  for  a  larger  computer  system  requires  only  a  few  changes  in 
the  program.)  In  order  to  represent  the  gravity  field  surface  suf- 
ficiently-well,  a  high  density  of  prediction  points  (grid  points)  has  to 
be  chosen.  If  not  more  than  100  data  are  found  in  the  prediction 
region,  the  algorithm  uses  all  the  data,  calculates  the  inverse  of 
the  corresponding  covariance  matrix  and  predicts  the  signals  tog  'ther 
with  its  rms-errors  at  all  grid  points.  The  situation  changes  if 
more  than  one  hundred  data  are  used.  The  algorithm  switches  over 
to  a  mode  which  can  best  be  described  as  moving  inverse  covariance 
prediction.  Its  principle  is  as  follows.  Assume  a  fairly  homogeneous 
data  distribution  and  a  grid  as  in  Fig.B2.1.  gravity  field  surface 
function  values  are  to  be  predicted  at  all  the  grid  points.  Assume 
furthermore  a  circular  region  R[  centered  at  the  grid  point  Pi  and 
another  circular  region  R^  centered  at  a  neighbouring  grid  point 
Pi+1  (see  Fig.  B2. 1).  All  data  within  Ri  are  used  for  the  prediction  of 
Si  (the  gravity  field  surface  function  value  at  the  grid  point  Pi),  all 
data  within  R^  are  used  for  the  prediction  of  Si+i,  have  a  common 
subset  Si|  i+i  which  is  the  intersection  of  Si  and  Si+1, 

Si,  i+1  =  Si  p  Si+i  . 

In  Fig.B2.l  the  subset  S^  2  consists  of  the  data  (black  dots)  within  the 
crosshatched  region. 


) 
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Figure  B2. 1:  Data  selection  in  moving  inverse  covariance  prediction 
mode 

Let  Ni  and  Ni+1  be  the  number  of  data  belonging  to  the  set  Si  and 
SyJ-1,  respectively,  and  let  Ni,  i+l  be  the  number  of  data  belonging 
to  the  common  subset  Si,  i+l.  (in  Figure  B2. 1,  Ni  =  28,  N2  =  30, 

Nit  2  =  27.)  The  corresponding  covariance  matrices  are  Ci  and  Ci+i, 
its  common  part  Ci,  i+l.  If  the  grid  is  dense  relative  to  the  data 
distribution  (a  necessary  requirement  for  surface  prediction),  then 
C i+l  will  differ  only  slightly  from  Ci;  with  other  words,  the  dif¬ 
ference  matrices  Ci  Ci,  i+l  and  Ci+l-Ci,  i+l  will  have  only  a  few  non¬ 
zero  elements.  (in  Fig.  B2. 1  ^  has  406,  C2  has  465,  Cl,  2  has  378; 


therefore  Ci-Cl(  2  has  only  28,  Cz~Cit2.  only  87  in  general  non-zero 
elements.  Here  we  have  considered  only  the  upper  triangular  part 
of  the  symmetric  covariance  matrix.  )  Consequentlv,  it  should  be 
possible  to  find  Cpfi  (the  inverse  covariance  matrix  corresponding 
to  point  Pi+i)  in  a  fast  and  simple  way,  if  is  known.  In  fact, 
the  problem  reduces,  in  principle,  to  a  two-fold  application  of  matrix 
inversion  by  block -partitioning. 

T 

Let  C|  be  partitioned  into  4  parts  Bi,  Bi  ,  and 

-1  T 

similarly  its  inverse  Cj  into  K^,  L^,  Lj  ,  M^, 


Ci  =  /Ci,i+l 


(2.1) 


There  it  follows  from  simple  matrix  algebra  (Faddejew  -  Faddejewa, 
1970,  p.  201  ff. ),  that 


<Di  -  ci,\+i  Bi)'1. 

'c"i!m  BiMi' 

C-.1  (I  -  b.lt  ) 

1,  l+l  11 


(2.2a) 


(2.  2b) 


(2.2c) 


with  I  denoting  the  unit  matrix.  In  the  same  way  as  above,  the  matrice 

C...  and  can  be  partitioned, 

i+l  l+l 


Ci+1  =/Ci,  i+l  Bi+1  V 


K...  L... 
i+l  itI 

T 

LT,.  M... 
i+l  i+l 


(2.3) 


Note  that  C.  ...  is  common  to  both,  C.  and  C...  and  so  is  the  inverse 
l,  i+l  i  i+l  ^ 

C  ,  If  this  common  inverse  is  known,  then  C.t,  can  be  found  in  the 
i,  i+l.  i+l 

way  described  above, 
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;w  *  <Dm  -  Bli  c'i!m  Bi«)'1 

(2.3a) 

...  =  -C"1  B  .  M .  . 

i+l  l,  i+l  i+l  i+l  , 

(2.  3b) 

i+l  =  C i,  i+l  ^  "  Bi+1  Li+1^  * 

(2.  3c) 

Therefore,  the  problem  of  calculating  C.^,  if  C.^is  known,  consists  in 
finding  C }  and  in  calculating  M.+j,  L._^,  K._^  following  the  equations 
above. 

-1  T 

Multiplying  equation  (2.2b)  by  M^  L.  ,  we  obtain 
L.M’!lT  =  -C'.1  .B.LT, 

ill  l,  l+l  l  i 

which  gives,  by  inserting  into  (2.2c),  the  inverse  of  the  covariance  matrix 

corresponding  to  the  common  data  set  S,  ,  expressed  by  the  known  four 

-1  1*  1+1 
submatrices  of  C .  , 

i 


-1  T 

C.  ...  =  K.  -  L.M.L.  . 

l,  i+l  l  ill 


(2.4) 


The  important  point  is  that  as  well  as  is  usually  very  small 

compared  with  C.  ;  (for  the  data  configuration  of  Figure  2.1,  Mj 
has  dimension  (1,  1),  M2  has  dimension  (3,3)  and  Cj,2  has  dimension 

(27,27).)  It  is  important  to  realize  that,  for  the  transition  from  C }  to 

_1  1 

Cf+i  ,  only  two  inversions  of  generally  very  small  matrices  with  dimensions 

of  M.  and  are  necessary.  (The  vector  and  matrix  multiplication  are 

relatively  inexpensive.  ) 

There  needs  still  one  problem  to  be  solved  which  was  tacitly  passed 
by:  the  matrix  C.*  has  to  be  re-ordered  according  to  a  permintation  vector 
whose  elements  point  to  the  data  within  the  prediction  circle.  This  pro¬ 
cedure  consists  mainly  of  logical  operations  and  is  very  fast;  it  is  accomp¬ 
lished  by  the  subroutines  BUBBLE  (which  makes  a  vector  bubble  sort, 
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generating  the  permintation  vector),  and  MATUMD  (which  does  the  actual 
matrix  reordering). 

The  principle  is  as  follows:  data  were  found  in  the  region  S., 

a  vector  V.  is  stored;  its  elements  point  to  the  data  in  S.,  the  inverse 

1  -1  1 
covariance  matrix  C.  is  given;  in  the  next  step  N.^  data  are  found  in 

the  region  S.+^,  a  vector  V.+^  is  generated  whose  elements  point  to  the 

data  in  S.,,.  N.  ...elements  are  common  to  V.  and  V..,,  but  the  cor- 

l+l  i,  l+l  l  l+l 

responding  sequential  locations  within  V.  and  V.+^  will,  in  general  differ. 

Therefore,  V.+^  passes  a  sorting  algorithm  which  "bubbles"  all  common 

N.  data  in  V  upward  such  that,  after  the  bubble  sort,  the  N.  ... 
i,  l+l  l+l  r  i,  l+l 

common  elements  occupy  the  first  N.  .+^  places  in  the  vector  V.^.  A 
pointer  vector  W.  is  generated  whose  elements  point  to  the  retained 
elements  in  V^. 


j 

Vi(j) 

vi+iO) 

before  sort 

Vi+i(j) 

after  sort 

Wi(j) 

1 

3 

7 

3 

1 

2 

17 

8 

8 

3 

3 

8 

79 

24 

4 

4 

24 

3 

45 

6 

5 

11 

41 

7 

2 

6 

45 

24 

79 

5 

7 

45 

41 

8 

13 

13 

Table  2.1  Example  of  bubble  sort  used  in  the  prediction  part  of  GSPP 
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This  vector  serves  as  a  permutation  vector  for  the  re-ordering  of 
the  inverse  covariance  matrix  C^"^. 

At  this  point  the  reader  might  ask  why  we  are  calculating 
the  inverse  covariance  matrix;  we  could  probably  apply  a 
similar  algorithm  for  the  calculation  of  the  solution  vector  of  the 
linear  system. 

C«  =  Cp  (2.  5) 

with  C. ..  covariance  matrix, 

CD. ..  cross-covariance  vector  (signal-data) 

r 

<y. . .  solution  vector 

Note  that  the  right  side  of  equation  (2.  5)  is  Cp  and  not  the  data 
vector  X;  this  can  be  done  since  we  estimate  for  each  solution  vector 
only  one  signal.  The  predicted  signal  is  then  given  by 

Sp  =  ft^X  (2.6  a) 

and  its  estimated  error  variance  by 

mp  =  CPP  "  aTcP  (2.6b) 

This  method  is  reportedly  two  to  three  times  faster  (Lachapelle,  1977) 
essentially  because  it  bypasses  the  matrix  inversion.  Here  the 
inversion  method  has  been  chosen  since  it  looked  more  transparent 
to  the  author;  the  solution  vector  method  should  bs  investigated. 

As  mentioned  before,  the  use  of  the  moving  inverse  prediction 
method  is  somewhat  restricted;  it  is  very  advantageous  for  the 
solution  of  very  large  problems  which  reduce  essentially  to  interpola¬ 
tions,  differentiations  and/or  downward  continuatious;  examples  are: 
determination  of  a  digital  geoid  based  on  altimeter  data,  solution  of 
the  Bjerhammar  problem,  prediction  of  mean  gravity  anomalies  from 
point  gravity  anomalies,  gravity  interpolation,  interpolation  of 
vertical  deflections,  determination  of  mean  gravity  anomalies  from 
gravity  and  g radiometer  data,  etc.  For  the  solution  of  problems 
which  involve  the  whole  data  vector  and  allow,  in  addition,  the 
estimation  of  parameters,  an  excellent  operational  system  is  available 
which  is  based  on  stepwise  least- squares  collocation  (Tscherhing, 

1974). 
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3.  Regression 

The  calculation  of  the  parameters  of  a  least- squares  regression 
polynomial,  in  GSPP,  is  based  on  the  following  premises:  the  data 
are  error-free  and  homogeneous  evaluation  functionals  of  a  surface; 
data  (and  surface)  are  defined  on  the  Euclidean  plane.  Then  a  least- 
squares  regression  polynomial  solution,  based  on  the  data  {x^3  .  {y.,3  , 
{f.3,  i=  1,  ....  I,  f  =  f(x  ,  y  ),  is  described  in  Section  2.3  of  Part  A; 
the  parameters  of  the  polynomial  are  given  by  equation  (2.3-2).  The 
design  matrix  $  =  {<$^3  =  £Lj  cpj }  ,  9j  ...  base-functions,  in  its 
explicite  form  ,  is  given  by 


*»♦ 

y». 

V2 

Xf> 

xtYi> 

yL 

n 

•  •  • »  Yl 

X2» 

Yz> 

X2 

x2» 

x2Y2» 

yf. 

n 

•  • •  t  Yz 

• 

• 

xr  yr 

X2 

*r 

xiyr 

* 

• 

n 

• •  •  >  Yj 

The  corresponding  normal  equation  matrix  (for  equal  weights)  is 
unstable  for  large  n,  GSPP  allows  the  degree  n  to  vary  between  0  and  5; 
the  degree  has  to  be  such  that  the  number  of  data  I  is  bigger  than  the 
number  of  parameters  J,  with  J  =  J  (n)  given  by 


J  * 


there  is  no  generalized  inverse  solution  allowed  in  GSPP.  If  n  has 

been  defined  such  that  1^  J  ,  the  highest  possible  degree  will  be  chosen. 

For  reasons  of  stability,  the  coordinates  (x.,  y.3  are  transformed 

i  i 

linearly  such  that  min(xl,  yp  =  0  and  nyix(x!  or  yp  =  1,  depending  on 
whether  the  range  of  x  or  y  is  bigger  .  Therefore,  also  the  para¬ 

meters  {a. }  ,  j  =  1,  ....  J,  refer  to  these  transformed  coordinates. 
3 
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This  has  to  be  taken  into  account  if  operations  of  any  kind  are  to  be 
applied  on  the  polynomial. 

The  actual  calculation  of  the  least- squares  regression  polynomial 
is  performed  in  the  subroutine  REGPOL.  If  the  domain  of  definition 
was  not  specified  by  the  user,  it  will  be  defined  in  GSPP  (if  REGPOL 
is  run  separated  from  GSPP,  these  values  have  tc  be  defined).  The 
output  is  the  vector  of  polynomial  coefficients,  normalized  a*  described 
above,  the  root  mean  square  and  average  absolute  approximation  error, 
the  individual  approximation  errors  (actual  function  values  minus  poly¬ 
nomial  derived  function  values)  if  requested,  and  a  matrix  of  polynomial 
derived  function  values  at  the  grid  points  of  a  user- specified  (or  GSPP 
selected)  regular  rectangular  grid.  These  grid  point  values,  in  turn, 
can  be  used  by  GSPP  for  profiling,  contouring  and  a  3-D  plot. 

In  the  following  we  give  some  examples  of  how  the  regression 
part  can  be  used. 

3.  1  Regression  polynomial  based  on  irregularly  distributed  data 

In  this  section  atypical  application  of  a  least- squares  regression 
is  shown:  there  is  given  a  set  of  error-free,  homogeneous,  irregularly 
distributed  data  defined  on  the  two-dimensional  Euclidean  plane.  A 
least- squares  regression  polynomial  of  a  certain  degree  has  to  be 
calculated  and  interpolated  at  the  grid  points  of  a  regular  rectangular 
grid.  In  addition,  the  individual  approximation  errors  (data  reproduction 
errors)  should  be  calculated.  In  the  following  program  the  polynomial 
degree  has  been  chosen  to  be  equal  to  4  which  is  too  high  relative  to 
the  number  of  data  which  was  selected  to  be  equal  to  14.  Therefore, 
the  program  changes  the  degree  to  3. 

In  the  sequel  the  program  plus  input/output  are  listed. 
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c 

c 

o 

0 

c 

0 

c 

c 

c 

c 

c 

c 

c 

t 

c 

(i 
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c 

c 

o 


o 

6C0Q 

0001 

C002 

>A)03 


AN  EXAMPLE  OF  A  LEAST-SQUARES  REGRESSION  P0LYN0MI  ALB  AS  ED  OR 
IRREGULARLY  DISTRIBUTED  DATA 
DIMENSION  A! 8,3, 1) ,  ZZ( 14) 

G0IEI0N  /:2AT/X(  14)  /YDAT/Y!  14)  /ZDAT/Z! 14)  /UGIP/IDUM1  ,MX,NY, 

#  IDUffiUO) ,  INI ,  IDUIISv  17)  ,  N  /UCFJVXLL,  YLL,  DX,  DY  /I0SPAR/I0UT 

A  STORES  THE  POLYNOMIAL  INTERPOL ATE D  FUNCTION  VALUES  AT  THE  GRID 
POINTS 

TJAf  IS  THE  NUMBER  OF  DATA 

IRR  ...  .CE.O  ...  DATA  ARE  IRREGULARLY  DISTRIBUTED 
U  . . .  DEGREE  07  POLYNOMIAL 

K( . ) .  Y(  .  >  ...  DATA  LOCATION  COORDINATES  (X  ACROSS (SOUTH  -  NORTH), 
Y  ALONG  -.VEST  -  EAST)) 

Z^ZvM.Y)  ...  FUNCTION  VALUES 

INI  ...  =1  OR  3  ...  ALSO  CALCULATION  OF  INDIVIDUAL  APPROXIMATION 
ERRORS:  IN  THIS  GASS  THE  ORIGINAL  FUNCTION  VALUES  WILL  BE  DES- 
TEOi'ED  AND  REPLACED  BY  THE  INDIVIDUAL  DATA  REPRODUCTION  ERRORS 
ICMAM,  N7MAX,  .121  LAX  ARE  THE  DIMENSIONS  OF  A  AS  DEFINED  IN  THIS 
PROGRAM 

MX,  NY  ...  NUMBER  OF  GRID  POINTS  IN  X  AND  Y  DIRECTION 
’  LZ.  Y!I  TY.  LIS.  UYrXAX 

xll/yllT. coordinates  of  the  grid’s  lower  left  comer  (corres¬ 
ponds  TO  A< 1 . 1 ,  1  > ) 

DX.  DY  .  .  .  GRID  DISTANCES  IN  X  AND  Y  DIRECTION 

I GUT  ...  OUTPUT  UNIT  FOR  MESSAGES  AND  RESULTS 

ICUIIl ,  IDUM2(G)  AND  IDUM3C  17)  ARE  DUMMY  VECTORS 

llDAT=  14 

IRR=0 

N=4 

nr:  - 1 
ilz:\::=q 

*  f  »— 

tixUrO 

mz:cx=  i 

HX=(. 

NY=C 

XLL=90. 

YLL= -130. 

1)11=4. 

DY2  •»  • 

nur=o 

READ  THE  DATA  (HERE  :  GENERATE  THE  DATA) 

DO  1  1*1, NDAT 
X(  I/=SIN!  I "2 . ) :S  10+ 100 
\\  I )  =  CCS (  Is.i3.Ki20- 100 
Z(  l!  =  CGSl  Ii.:0.U):.-!SIII(  1*0.9) 

2?J  I)  =Z(  I) 

CALL  THE  SUBROUTINE  REGPOL 

CALL  REGPOL l  A.  MXTIAX.  NY? LAX,  MZMAX,  NDAT,  I HR) 

PRINT  TIC  DATA  AND  RESULTS 
WRITS! IOUT.GOOO, 

DO  3  1=1, NDAT 

WRITE( I0UT,  6001)  I,X< I) ,Y( I) ,ZZ( I) ,Z( I) 

WRITE( I OUT, 0002) 

DO  3  1=1, MX 
Il=HX-l+l 

WRITE! I0UT, 6003)  CA(  1 1 , J, 1) ,J=1,NY) 

FORMAT!  1  1 4X»  •  I  X(I)  YC  I)  Z(I)  DZ(I)',//> 

FORILAT(  III  ,411. 12, 2C  UX,  2F3 . 2) ) 

FORMAT*  1.-0,  IX) 

FORMAT!  IliO,4X,3F8.2) 

STOP 

F.ND 


Program  corresponding  to  Table  3.  1.  1 


-130- 


I 


1 

2 
r\ 
O 

4 

5 
0 
7 

o 

9 

10 

11 

12 

10 
14 


0.73  0.57 
0.10  -0.04 
0.15  0.01 
0.47  0.23 
0.63  0.36 
0.22  -0.17 


X(I>  Y(I) 


100.09  -119.80 
92.43  -80.80 

97.21  -113.22 
109.89  -83.12 

94.56  -115.19 
94*.  63  —86.79 

109.91  -110.95 
97.12  -91.52 

92.49  -105.34 
109.13  -96.91 

99.91  -100.27 
90.94  -102.56 
107.63  -94.67 

102.71  -103.00 


0.34 

0. 10 

-0.22 

-0.40 

-0.  17 

-0.33 

0.07 

-0.11 

O.OS 

-0.  15 

-0.56 

-0.39 

-0.09 

-0.52 

-0.42 

-0.20 

-0.30 

-1.12 


Z(I)  DZ(I) 


0.69 

0.  18 

0.53 

-0.  13 

O.OS 

-0.43 

0.  13 

-0.03 

0.73 

0.40 

0.77 

0.33 

-0.02 

-0.33 

-0.52 

-0.49 

-0.20 

0.08 

0.  12 

0 . 43 

-0.S2 

-0.08 

-0 . 94 

-0.20 

-0.74 

-0.27 

0.03 

0.00 

-0. 19 

-0.  15 

0.08 

-0.54 

-0.40 

-0.06 

-0.09 

-0.20 

0.21 

1 

o 

»-* 

0.04 

0.48 

-0.30 

-0.  10 

0.33 

-1.19 

-1.05 

-0.66 

y 


Table  3.  1.  1:  Third  degree  least- squares  regression  polynomial 

based  on  irregularly  distributed  data  x;  y;  z  =  z(x,  y) 


3. 2  Regression  polynomial  based  on  regularly  distributed  data 

If  data  are  distributed  on  a  regular  rectangular  grid,  and  stored 
on  an  array  corresponding  to  this  grid,  then  a  slightly  modified  version 
of  the  above  listed  program  is  necessary  in  order  to  obtain  the  regression 
field.  The  following  example  may  serve  as  an  illustration. 

Both  regressions  described  above  can  also  be  obtained  by  .ailing 
GSPP;  the  regression  surface  can  be  contoured,  profiled  or  plotted 
as  a  3-D  view.  These  integrated  applications  will  be  described  in 
PART  C. 


G  REGRESSION'  POLYNOMIAL  BASED  ON  REGULARLY  DISTRIBUTED  DATA 
DIMENSION  A(3, 8, 2> 

COMMON  /UC  I?/  IDUMl ,  MX,  NY,  IDUM2(6)  ,  IW1, IDUK3C  17)  ,N  /UCRP/XLL,  YLL,DX 
tt  , DY  / lOCPAR/ I OUT 

NDAT=0 
IRR--1 
K=5 
IWi=l 
1221411=8 
NYKAX=3 
ME?' AT '=2 
121=6 

I  IVs  8 
XLL=90. 

YLL= -120. 

BX=4. 

DY=  !5 . 

ICUT=6 

G  READ  THE  DATA;  HERE:  GENERATE  THE  DATA 
BO  1  1=1,  MX 
BO  1  J= 1 , NY 

1  A(  I ,  J ,  1 )  =  C  S  INC  1*2 . )  *COS(  J*3 . )  )  *  10 
C  PRINT  THE  DATA 

BO  4  1=1, MX 

I I  =  121- 1 -M 

4  WRITE! I OUT, 6003)  (  A( 1 1 , J, 1) , J= 1 ,NY) 

G  CALL  THE  SUBROUTINE  REGPOL 

CALL  REG?OL(  A,  IEGIAX,  NYJ1AX,  MZMAX,  NDAT,  IRR) 

C  PRINT  THE  GRID  VALUES  0?  THE  REGRESSION  POLYNOMIALS  AND  THE  INDI- 
C  VI  DUAL  APPROXI HAT I  ON  ERRORS 

DO  2  1=1, MX 
I1  =  HX-I+1 

WRITE! I OUT, 6003)  (A( II, J, 1) ,J=1,NY) 

2  CONTINUE 

BO  3  1=1, MX 
1 1  =  121-1+1 

WRITE!  I OUT, 6003)  ( A( 1 1 , J, 2) , J= 1 ,NY) 

3  CONTINUE 

6003  FORMAT! 1K0, 4X.8F8.2) 

STOP 

END 


Program,  corresponding  to  Table  3.2.  1 
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5.31 

-5.  15 

4.89 

-4.53 

5.39 

-5.22 

4.96 

-4.59 

-9.79 

9.50 

-9.01 

8.35 

2.77 

-2.63 

2.55 

-2.36 

7.49 

-7.27 

6.90 

-6.39 

-9.00 

3.70 

-8.23 

7.67 

Regularly  distributed  data 


4.03 

•-3.34 

2.94 

-2.28 

4. 13 

-3.59 

2.98 

-2.31 

-7.52 

6.53 

-5.42 

4.20 

2. 12 

-1.85 

1.53 

-1.19 

5.75 

-5.00 

4. 15 

-3.21 

-6.91 

6.00 

-4.98 

3.36 

AVERAGE  ABSOLUTE  POLYNOMIAL  APPROXIMATION  ERROR  . . 
RII3  POLYNOMIAL  APPROXIMATION  ERROR  ...  5.483 


4. 804 


4.06 

-0.55 

1.71 

-0.48 

•1.53 

-2.  12 

0.25 

0.01 

2.  19 

2.39 

•4.32 

-0.83 

-1.02 

-0.40 

-0.08 

0.64 

-1.25 

-0.44 

0.  16 

0.01 

1.75 

0.37 

0.49 

0.24 

0.00 

-0.05 

0.36 

0 . 38 

0.02 

0.36 

■0.21 

-0.03 

•0.97 

-1.31 

0.30 

-0.07 

-0.39 

-1. 11 

-0.09 

-1.54 

0.68 

0.26 

0.66 

0.86 

-1.14 

-1.  12 

1.04 

1.35 

4th  degree  regression  polynomial  grid  point  values 


x 

A 

— >y 


1.25 

-4.60 

3.67 

-4 . 74 

8.26 

11.62 

2.52 

-2.69 

5.30 

-9.66 

•4.  18 

9.61 

5.91 

-4.  13 

5.04 

-5.23 

-7.76 

8.79 

2.38 

-2.37 

5.  15 

-6.75 

-8.73 

7.44 

4.07 

-3.49 

3.27 

-4.  17 

-7.54 

6.  17 

2.33 

-1  .82 

6.72 

-3.49 

-6.60 

6.08 

3.33 

-1.17 

3.07 

-0.76 

-6.  iO 

3.94 

0.87 

-2. 04 

3.28 

-2.09 

-6.02 

2.31 

Residuals 

Table  3.  2.  1:  Least- squares  regression  polynomial  based  on  regularly 


distributed  data. 
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4.  Smooth  surface  representations 

One  essential  part  of  GSPP  is  a  set  of  subroutines  which  is 
needed  for  a  smooth  representation  of  a  surface  by  a  bicubic  spline 
function  together  with  a  surface  interpolation/differentiation  algorithm. 

4.1  Calculation  of  spline  defining  values 


This  section  gives  a  practical  example  of  how  to  obtain  the  spline 
defining  values  by  the  subroutine  BISP,  isolated  from  GSPP.  The  neces¬ 
sary  background  and  formulas  are  contained  in  Section  2.6.2.  of  PART  A. 

Let  us  assume  that  function  values  are  given  at  the  grid  points 
of  a  regular  rectangular  grid  (m  =  l,  ...,  M;  n  =  l,  . . . ,  N).  Then  the 
bicubic  spline  defining  values  is  the  set  jfj^,  pmn,  Imn*  rmuf  »  m=l» 
...,  M;  n=l,  ...»  N  (see  Section  2.6.2)  of  PART  A.  The  function  values 
at  the  grid  points  |fmn^  are  assumed  to  be  known,  the  derivatives 
{  Pmn’  ^mn*  rmn  I  have  to  be  determined.  Under  side  conditions  as 
explained  in  Section  2.6.2,  of  PART  A  this  set  of  values  defines  the  bi¬ 
cubic  spline  surface  uniquely. 

The  following  program  can  be  used  to  calculate  these  values  by 
using  BISP. 


C  AN  EXAMPLE  OF  A  BICUBIC  SPLINE  DEFINING  VALUE  CALCULATION 
DIMENSION  AC  10. 10,4) 

C  AC ....  1)  ...  STORES  TIIE  REGULARLY  DISTRIBUTED  DATA 

C  AC . , . ,2)  . . .  1.  DERIVATIVES  IN  X-DIRECTION  CACROSS,  SOUTH  -  NORTH) 

C  AC.,., G)  ...  1.  DERIVATIVES  IN  Y-DIRECTTON  (ALONG,  NEST  -  EAST) 

C  AC.,,, 4)  ...  2.  HIKED  XY-DERI VATIVES 

C  NX,  NY  ...  ACTUAL  DIMENSION  OF  THE  GRID 
C  OUTPUT  UNIT  =  <5 
HX=6 
NY=0 
I0UT=6 

C  READ  THE  DATAi  HERS:  GENERATE  THE  DATA 
DO  1  1=1,  IK 
DO  1  J= 1 , NY 

1  AC I,J, 1)=<SINC0.G*I)+C03(0.6*J) )*10 

C  CALCULATE  'HIE  DEFININC  VALUES 

CALL  BISPCA, 10, lO.IK.NY) 

C  PRINT  DATA  AND  RESULTS 

DO  G  K=  1 ,4 
VRITEC I OUT, 0900) 

DO  2  1=1,  IK 
1 1=HX-I+ 1 

2  WRITE!  IOUT.OOOl)  C  AC  1 1 > J ,K) , J= 1 , NY) 

3  CONTINUE 

<5000  F0RIIATC 1  HO ,  1 X,  /V) 

600 1  FORMAT!  1  IiO ,  4X,  5F 1 0 . 4) 

STOP 

END 


Program  to  Table  4.  1.  1. 
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Table  4 

.  1. 1 

Bicubic 

spline  defining  values 

4 

: 

9.6646 

•  5.0348 

-0.8608 

-5.9627 

-8.4887 

3  • 

14.2381 

9.6083 

3.7127 

-1.3392 

-3.9152 

Function 

values  at  the 

17.S46& 

12.7166 

6.S210 

1.7190 

-O.SO70 

grid  points  (data) 

X 

18.2283 

13.5985 

7.7029 

2.6010 

0.0750 

/ 

V 

16.6681 

■•o  (•scrn 
Iw  •  WOuu 

6.1427 

1.0408 

- 1 . 4352 

- >y 

1 

13.6476 

8.4178 

2.5222 

-2.5797 

-5 . 1057 

t 

! 

* 

-4.3432 

-4.8432 

-4.8432 

-4.3432 

-4.3432 

i 

-4.0341 

-4.0341 

-4.0341 

-4.0341 

-4.C041 

Dx-derivatives 

1 

-2. 0637 

-2.0657 

-2.0657 

-2.0657 

-2.0657 

1 

0.3262 

0.3262 

0 . 3262 

0.3262 

0.3262 

2.7953 

2.7953 

2.7958 

2.7958 

2.7953 

4.0323 

4.0328 

4,0323 

4.0328 

4.0028 

1 

-4.2800 

-5.3293 

-5.9790 

-3.7474 

-1.9153 

f 

-4.2300 

-5 . 3293 

-5.9790 

-3.7474 

-1.9153 

-4.2E00 

-5.3293 

-5.9790 

-3,7474 

-1.9153 

Dy-derivatives 

-4. 2800 

-5.3293 

-5,9790 

-3.7474 

-1.9153 

-4.2800 

-5.3293 

-5.9790 

-0.7474 

-1.9153 

-4.2300 

-5 . 3293 

-5.9790 

-3.7474 

-1.9153 

) 

-0.0000 

0.0000 

0.0000 

-0.0000 

-0.0000 

-0.0000 

-0.0000 

0.0000 

-0.0000 

0.0000 

Dvv-denvatives 

0.0000 

0.0000 

-0.0000 

0.0000 

-0.0000 

I 

0 • CG00 

0.0000 

Q.CQC0 

-0.0000 

-0.0000 

1 

I 

1 

-0.0000 

-o.ccoo 

-0.0000 

-0.0000 

0.0000 

1 

» 

1  ■ 

0,0000 

o.ccoo 

-O.COOO 

0.0000 

0.0000 

i 

i  i 

1 


-135- 


The  reader  will  probably  realize  that  no  grid  distances  have 
been  defined.  This  is  not  necessary  because  the  defining  values 
returned  from  BISP  are  normalized  ones:  they  refer  to  a  square 
grid  of  grid  distance  equal  to  1.  The  reasons  for  the  normalized 
calculation  are  simplicity  and  stability;  there  is  no  restriction  of 
generality  involved:  it  can  easily  be  shown  that  the  defining  values 
referring  to  a  non -normalized  grid  can  simply  be  obtained  by  dividing 
the  x-aerivatives  by  the  grid  distance  in  x-direction,  the 

y-derivatives  Iq^m)  by  the  grid  distance  in  y-direction  and  the  second 
mixed  derivatives  Crmn3  by  the  product  of  x-  and  y-  grid  distances. 

BISP  is  designed  for  a  maximum  grid  size  of  300  x  300.  The 
CPU-time  needed  increases  only  linearly  with  the  number  of  grid 
points  involved,  a  good  rule  of  thumb  is:  number  of  grid  points  # 

1.2  •  10'4  seconds;  this  number  refers  to  a  AMDAHL  470  V/6  -II 
computer.  The  following  Table  lists  CPU-  estimates  for  a  couple  of 
grid  sizes  n  (square  grid).  The  estimates  refer  to  a  AMDAHL  470 
V/6-II  computer: 


n 

CPU-time  (sec) 

5 

0.005 

10 

0,012 

25 

0.074 

50 

0.288 

100 

162 

Table  4.  1.  2 


4.  2  Smooth  surface  interpolation /differentiation 

How  a  bicubic  spline  is  interpolated  and  differentiated, 
is  demonstrated  in  subsection  2.6.3;  the  formulas  (2. 6 -7  a,  b)  are 
optimal  in  terms  of  computer  time  (PART  A). 

The  interpolation/differentiation  of  a  spline  is  a  3- stage  process: 
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In  the  first  step  the  grid  element,  corresponding  to  the 
coordinates  of  the  calculation  point,  has  to  be  found:  with  (x*,  yj)  the 
coordinates  of  the  lower  left  grid  point  and  (h^.,  h  )  the  grid  distances 
in  x-  and  y-  direction,  the  grid  element  indices  are,  for  a  point  (x,  y), 
simply  given  by 

m  =  int((x-Xj)/h^)  +  1,  n  =  int  ((y-yj/hy)  +  1, 

(int.  ..  integer),  if  the  calculation  point  (x,  y)  is  located  within  the 
domain  of  definition  of  the  bicubic  spline. 

The  second  step  consists  in  the  calculation  of  the  16  coefficients 
of  the  bicubic  element.  A  straightforward  way  would  be  a  calculation 
using  equations  (2.6 -5a,  b);  which  gives  the  4x4  matrix  of  coefficients 
as  a  product  of  three  4x4  matrices  (see  PART  A). 

A  =  HT(hx)FH(hy)  . 

This  operation  involves  128  multiplications  and  128  additions  and  uses, 
even  in  the  normalized  version  (h^  =  hy.  =  l),  700  micro- seconds  to 

calculate  all  16  coefficients  (the  matrix  A).  Since  A  has  to  be  cal¬ 
culated,  in  general,  for  each  calculation  point,  it  is  very  important  to 
optimize  this  algorithm  in  terms  of  CPU- time.  After  many  trials  I 
found  a  very  fast  and  probably  optimal  solution  which  involves  no 
multiplication  and  only  58  additions;  the  time  elapsed  for  the  calculation 
of  all  16  coefficients  using  this  fast  algorithm  is  as  little  as  58  micro¬ 
seconds  which  corresponds  to  a  12-fold  gain  in  calculation  speed. 

The  following  Table  lists  CPU-estimates  for  the  calculation  of 
all  16  coefficients  for  all  elements  of  varying  size  square  grids 
(nxn  elements).  The  estimates  refer  to  a  AMDAHL  470  V/6-II  computer: 


n 

CPU-time  (sec) 

5 

0.0015 

10 

0.0058 

25 

0.0363 

50 

0. 1450 

100 

0.5800 

TABLE  4.2.  1 


-137- 


The  subroutine  BILDE  is  responsible  for  the  calculation  of  the 
coefficients. 

After  these  two  steps  the  actual  interpolation/differentiation  can 
be  performed.  Since  all  calculations  in  the  bicubic  spline  algorithms 
refer  to  a  normalized  grid  of  unit  grid  distances,  the  relative 
coordinates  of  the  calculation  point  are  also  to  be  referred  to  this 
unit  grid  distance;  for  a  point  (x,  y)  they  are  simply  given  by 

x  =  (x-xj)/!^  -  (m-i),  y  =  (y-yt)/hy  -  (n-t) 

with  (m,  n)  element  indice"  as  defined  above.  The  interpolated/diff  erentialed 
bicubic  spline  at  a  point  (x,  y)  (within  the  domain  of  definition)  can  then 
be  obtained  by  the  set  of  formulas  (2.6-7a,  b)  of  PART  A  after  division 
by  the  appropriate  grid  distances.  The  calculation  itself  is  performed 
in  the  subroutine  BSFC.  The  subroutine  is  designed  such  that  it  can 
provide  all  spline  derivatives 


dx 


a.,  a, 

‘ay  “ 


os  a*  +  a2  5  4,  ttj,  a2  s  2 


(BSFC  is  not  designed  for  third  order  derivatives  with  respect  to  \ 
and/or  y. ) 

The  following  Table  gives  a  listing  of  CPU-time  estimates  for 
the  interpolation /differentiation  part  and  for  the  total  CPU-time  used 
(index  finding,  calculation  of  parameters,  interpolation/differentiation). 


The  following  Figure  4.2.  1  with  the  corresponding  Table  4.2.3 
shows  the  function  values  of  25  regularly  distributed  data,  the  corres 
ponding  spline  surface,  and  lists  interpolated/differentiated  values. 

The  program  below  has  been  used  to  generate  the  output  given  in 
Table  4.2.3. 
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A 

0 
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6000 

6001 

6C02 

6000 


A  BICUBIC  SPLINE  I NTERPOLATION/D I FFERENTI ATI  ON  EXAMPLE 

A< . )  ...  ARRAY  STORING  THE  BICUBIC  SPLINE  DEFINING  VALUES 

HXMAX,  IIYMAX  ...  MAXIMUM  1.  AND  2.  DIMENSION  OF  A 
IK,  UY  ...  ACTUAL  1.  AND  2.  DIMENSION  OF  A 

SX.  DY  ...  GRID  DISTANCES  IN  X  AND  Y  -  DIRECTION  (ACROSS,  ALONG) 

XL,  XU  ...  LONER  AND  UPPER  X  -  BOUNDS  OF  THE  DOMAIN 

YL,  YU  ...  LONER  AND  UPPER  Y  -  BOUNDS  OF  THE  DOMAIN 

NC  ...  NUMBER  OF  CALCULATION  POINTS  PROCESSED 

X,  Y  . .  COORDINATES  OF  THE  CALCULATION  POINT  (ACROSS,  ALONG) 

XNR,  YNR  . . .  NORMALIZED  RELATIVE  COORDINATES  FOR  BSFC 
IX,  IY  ...  NUMBER  OF  DERIVATIVES  IN  X  AND  Y  -  DIRECTION 
DIMENSION  A( 5,5, 4) 

COMMON  /C0M2/XNR,  YNR,  DX,  DY 

H3EIAX=5 

NYHAX=5 

I!X=5 

NY=U 

DX=2.o 

DY=4. 

(IX. XU;  YL, YU)  ARE  THE  DOMAIN  LIMITS 
XL=1. 

YL= - 1 . 

XU=  XL*  ( MX- 1 )  SDX 
YU-  YLt  ( NY- 1 )  »DY 

PRINT  DOMAIN  AND  GRID  INFORMATION 
NRITEt 6 , 6000)  XL, XU,  YL,  YU , DX,  DY,  MX, NY 

READ  THE  BICUBIC  SPLINE  DEFINING  VALUES;  HERE  THEY  ARE  GENERATED 
DO  1  1=1,  MX 
DO  ■  J= 1 ,NY 

A(  I ,  J .  1)  =  (SIN(  I*1.)+C0S(J*2.)>*10 
CALL  3 IS?(  A,  MXKAX,  NYMAX. MX, NY) 

READ  THE  NUMBER  OF  CALCULATION  POINTS 
RE  AO  (5,  so  NC 
NR!TE( 0,6002) 

M=  0 
N=0 

BO  2  1=1, NC 

READ  ALL  CALCULATION  POINT  COORDINATES  AND  THE  NUMBER  OF  DIFFEREN¬ 
TIATIONS  IN  X( ACROSS)  AND  Y( ALONG)  DIRECTION 
READ; 3, s)  X.Y.  IX,  IY 

IF  A  CALCULATION  POINT  HAPPENS  TO  BE  OUTSIDE  THE  DOMAIN,  PRINT  A 
MESSAGE  AND  GO  TO  THE  NEXT  POINT 

IFCX.GE.XL. AHD.X. LT.XU.AND. Y.GE. YL. AND. Y.LT. YU)  GOTO  3 
NRITE(O.OOOO)  I, X.Y, IX, IY 
GOTO  2 
CONTINUE 

CALCULATE  TEE  EISMENT  INDICES  AND  THE  NORMALIZED  COORDINATES 
>111=  IX- XL) /DX 
MII=  INTI 'GO 

:e!f.=xn-mn 


MN=HN+1 
YN= ( Y-YL) /DY 
NN=  !NT(  YIO 
YNR=YN-NN 
NN=NN+ 1 

CHECK  IF  THE  INDICES  HAVE  CHANGED  RELATIVE  TO  THE  LAST  CALCULATION 
POINT;  IF  NOT.  THE  POLYNOMIAL  COEFFICIENTS  NEED  NOT  BE  CALCULATED 
AGAiH 

IFCM.EQ.HH.AND.N.EQ.NN)  GOTO  4 
11=121 
51=  UN 

CALCULATE  THE  COEFFICIENTS  OF  THE  BICUBIC  POLYNOMIAL 
CALL  B ILDE( A, MXMAX, NYMAX, M.N) 

CONTINUE 

I NTERF OL AT I ON/D 1 FFERENT I AT I ON 
F=ESFC(  JX,  IY) 

PRINT  THE  RESULTS 
NRITE(6,6001)  I, X.Y,  IX,  IY, 7 
CONTINUE 

FORMAT! ISO, AX, i3,2X,2F10.2,3X,2I3,2X, •  OUTSIDE  DOMAIN’) 

FORMAT!  I  HO,  AX,  I5,2X,2F10.2,3X,2I3,2X,F10.2) 

FORMAT! ISO, AX, ’  IX  Y  IX  IY  F’,// 

*) 


FORMAT! 1H1.AX, ’ 
•#.’ ,?10.2,5X,  ’XU  . 
a/, OX,  ’DX  . . .  ’  ,710 


DOMAIN  LIMITS  AND  GRID  PARAMETERS’ ,//,9X, ’XL  . 
’ ,F10.2,/,9X, ’ YL  . . . ’ .F10.2.SX, ’YU  ...’,  F10.2 
,5X,’DY  . . . * ,F10.2,/,9X, ’MX  . . . ’ , 1 10, 3X. ’NY  .. 


a’ ,  no,///) 

STOP 


END 


Program  corresponding  to  Table  4.2,3 


only  inter- 

total  CPU-time 

0(| 

0% 

pol  /diff. 

(micro sec. ) 

0 

0 

18 

76 

1 

0 

19 

77 

0 

1 

19 

77 

2 

0 

18 

76 

1 

1 

21 

79 

0 

2 

28 

76 

2 

1 

20 

78 

1 

2 

20 

78 

2 

2 

19 

77 

Table  4.  2.  2 

CPU  -  time  estimates  for  2-D  spline  interpolation /differentiation 


DOMAIN  LIMITS  AND  GRID  PARAMETERS 


XL  *  *  , 

1. 00 

xu  ... 

11.00 

YL  ... 

-1.00 

YU  .. . 

15.00 

DX.  . . . 

2.50 

DY  , .  . 

4.00 

NX  . . . 

5 

NY  ... 

5 

I 

X 

Y 

IX 

IY 

F 

1 

3.27 

7.38 

0 

0 

19.  12 

2 

8.53 

2.11 

l 

0 

-2.20 

3 

1,16 

10.24 

0 

1 

-4.04 

4 

13.  15 

7.55 

0 

OUTSIDE  DOMAIN 

5 

1.10 

14. 12 

I 

1 

0.00 

6 

9.99 

.2.53 

0 

o 

2.31 

7 

6.34 

12.2! 

2 

1 

0.00 

8 

7.51 

3.77 

1 

2 

-0.00 

9 

-7.22 

10.30 

2 

o 

OUTSIDE  DOMAIN 

Table  4.2,3:  Bicubic  spline  interpolation/differentiation 
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Figure  4.  2 


TfiTo  i 

8.0 


1:  Data  distribution  and  corresponding  bicubic  spline 
surface  contour  plot. 


1 1  v  * ' 


j.v.1  W'J 


.'.id 


Cl-'' 


-141- 


4.  3  Data  approximation  errors 

Whenever  the  function  values  of  the  surface  are  predicted  at 
the  grid  points  of  a  regular  rectangular  grid,  and  the  surface  is 
represented  by  a  bicubic  spline,  interpolating  these  grid  values,  the 
irregularly  distributed  homogenous  data  are,  in  general,  not  repro¬ 
duced.  It  is  obvious  that  the  data  reproduction  errors  (the  difference 
between  the  actual  data  f j  and  the  corresponding  surface  value  f  ^) 
decrease  with  decreasing  grid  distance;  if  the  error  norm  is  defined 
in  terms  of  the  maximum  error, 

j | e 1 1  :=  max  |fi  ••  fi|, 

Vi 

the  error  norm  goes  to  zero.  This,  however,  does  not  necessarily 
mean  that  the  spline  representation  of  the  whole  surface  is  getting 
better  with  decreasing  grid  distance.  Nevertheless,  the  above  error 
norm  provides  a  first  estimate  of  how  well  the  data  are  reproduced. 

The  subroutine  REPRO  is  designed  for  the  calculation  of  the 
individual  data  reproduction  errors,  the  average  absolute,  and  the  root 
mean  square  (RMS)  approximation  error.  When.ever  the  individual 
errors  are  calculated,  the  original  function  values  will  be  destroyed 
and  replaced  by  the  corresponding  approximation  (reproduction)  errors. 
If  a  data  point  happens  to  be  located  outside  the  surface  domain,  it 
will  not  be  considered  for  the  data  reproduction  error  calculation  and 
a  value  of  99999.99  will  be  assigned  to  its  error.  The  average  and 
RMS  errors  are  edited  on  the  lineprinter  (or  any  other  assigned  output 
device).  The  following  example  may  illustrate  an  isolated  use  of  the 
subroutine  REPRO.  Further  information  can  be  found  in  the  comment 
statements  to  REPRO* 


C  TEST  OF  THE  DATA  REPRODUCTION  ERROR  CALCULATION  PROGRAM  ’REPRO* 
DIMENSION  A  (25 1 25  » 4 ) »  ZZ(IOC) 

COMMON  /XOAT/XUOO)  /YDAT/Y(100)  /ZOAT/ZdOO)  /UCIP/IUPt 30) 
*/UCRP'0P(30>  / IOSPAR^NOUT 
DATA  HXMAX*NYMAX»MZMAX/2*25»4/ 

C  INPUT  UNIT  ...  INPUT*  OUTPUT  UNIT  ...  NOUT 

1NPUT=5 
NOUT-6 

C  READ  THE  DATA 

C  READ  THE  DATA*  ERROR  INDICATES  END  OF  DATA  Self  COPY  Z(.l 

C  ONTO  ZZ(.) 

N  =  1 

18  READ(INPUT»**£RR=19)  X(N)»  Y ( N) ♦  Z(N) 

ZZ(N)=Z(N) 

N-N+l 
bOTO  18 

19  N0AT=N-1 
WRlTEtNOUT *5001 >  NDAT 

5001  FORMAT( 1 HO* 4X ♦ ’NUMBER  OF  OATA  FOUND  IN  OATA  SET  ...  ’»I5»/> 

C  READ  AND  ECHO  SURFACE  ARRAY  PARAMETERS 

READ ( 5 » *)  IUP( 2 )  *  IUP(3)»  ( UP ( I ) *1  =  1 1 4) 

WRIT£(NOUT» 5002)  ( UP ( I ) *  1  =  1  *  4) »  IUP(2)»  IUP(3) 

5C02  FORMAT ( 1HQ *4X * ’SURFACE  ARRAY  PARAMETERS  :*./*5X» 

* ’COORDINATES  OF  LOWER  LEFT  SURFACE  ARRAY  POINT  1  X  =  »*F10.2* 

*•  Y  =  ’*F10.2*/*5X*’GRI0  DISTANCES  IN  X-ANQ  Y-D1RECT  ION  I*. 
*2F10.4*/»5X* ’NUMBER  OF  GRID  POINTS  IN  X-AND  Y-OIRECT ION  :*» 
*215*") 

C  READ  THE  SURFACE  ARRAY  A(.*.»l)  AND  ECHO  IT 

REAOC 5 » »  1  ( <A<  I#J*1>  *J  =  1» IUP131  >» I=IUP(2>  *1 »-l) 

WRITE (NOUT# 6000) 

60U0  FORMAT  11 HO *4X#*FU NOTION  VALUES  AT  ThE  GRID  POINTS  :•»//> 

CALL  ORUCK(.'\»MXHAX»NYHAX*MZMAX»  IUP(2>  #IUP<3> ♦ 1 »UP< 1 ) »UP( 2 ) *UP( 3 ) » 
•UP (  4  >  *  36  *8  *0  *0  *0 ) 

C  CALCULATE  THE  BICUBIC  SPLINE  REPRESENTATION 
CALL  BISP(A*HXMAX»NYNAX*IUP(2) »  IUP(3>) 

C  CALCULATE  THE  OATA  REPRODUCTION  ERRORS 
CALL  REPROU»MXHAX*NYMAX*NOAT»  3) 

C  PRINT  THE  DATA  AND  THE  INDIVIDUAL  ERRORS 
URITEtNOUT  »7u00 ) 

7000  F OR MAT(1 HO »4X** SPLINE  OATA  REPRODUCTION  ERRORS  I’*/* 

•5  X  *  *  I  X 1 1>  Y ( I )  Z(I)  ERRtl)**') 

DO  7001  1 1=1  *  NDAT 

7001  URITElwOUT  *7002 )  II*  X(I1>*  Y(II>»  ZZ(II>»  Z(II) 

7002  FQRMATUH  *4 k ♦  1 5 *  4F 10 . 2  ) 

STOP 
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NUHBtR  Of  OAT*  POUNO  I  ft  OAT*  SET  ...  17 

SURFACt  ARRAY  PARAMET EKS  : 

LOORulNATSS  OF  LOWER  LEFT  SURFACE  ARRAY  POINT  ;  X  r  *7.67  Y  =  14.00 


I>R10  01STANCES 

IN  X-ANO  Y- 

DIRECTION 

:  .0*17 

.0667 

NUMBLR 

OF 

bRIO 

POINTS  IN  X 

-ANO  Y-01RECT ION  : 

a  a 

rUNCT ION 

VALUES 

AT  THE  GRID  POINTS  : 

long. 

1* 

.0 

1*  *.0 

1*  8.0 

1*  12.0 

1*  16.0 

1*  20.0 

1*  2  *  .0 

1*  28.0 

lat. 

*7 

57.5 

-11.3V 

-12.17 

-7.12 

-12.16 

-11 .53 

-8.71 

-10.20 

-10.15 

*7 

55.0 

-10.25 

-*.21 

■'9.0* 

-11.25 

-11.53 

-9.78 

-10.73 

-10.96 

*7 

52.5 

-13.25 

-13.83 

-1*  .62 

-1* .63 

-11.81 

-12.33 

-12.01 

-13.79 

*7 

50.0 

-13.26 

-13.23 

-11.79 

-11.55 

-15.02 

-10.69 

-12.37 

-12.05 

*7 

*7.5 

-13.51 

-15.9* 

-18.65 

-28.03 

-20.28 

-20.53 

-20.01 

-18.73 

*7 

*5.0 

-18.30 

-19.75 

-20.72 

-19.62 

-28. *9 

-18.33 

-15.81 

-14.29 

*7 

*2.5 

-17.15 

-17.06 

-19.26 

-2*. *5 

-19.** 

-19. *5 

-19.78 

-19.5* 

*7 

*0.0 

-18.62 

-18.79 

-19.31 

-2*. *8 

-22.72 

-18.6* 

-16.89 

-17.4* 

AVERAGE  APPROXIMATION  ERROR  ...  .723 

RMS  APPROXIMATION  ERROR  ...  .8*8 


SPLINE  DATA  REPRODUCTION  ERRORS  1 


1 

XII) 

Y  A I  > 

2(1) 

ERR(I) 

1 

*7.88 

1  *.  26 

-11.35 

.*6 

2 

*7.8* 

1*  .26 

-13.*6 

.93 

3 

*7.77 

1*.17 

-22.93 

-.03 

* 

*7.72 

1*  .33 

-20.6* 

-1.63 

5 

*7.93 

1*.25 

-12.36 

-.53 

6 

*7.96 

1*.01 

-18.02 

lOUOOO .00 

7 

*7.70 

1  *.  16 

-21.37 

.93 

8 

*8.15 

1*.*3 

-11.33 

100000.00 

9 

*8.19 

1  *.36 

-1*.02 

iouoqo. oo 

10 

*7.75 

1*  .3* 

-17.0* 

.*6 

11 

*7.79 

1  *.23 

-2*. 51 

1.26 

12 

*7.93 

1*  .07 

-9.38 

-.20 

13 

*7.92 

1*.  52 

-10.70 

100000. 00 

1* 

*8.16 

1**53 

-10.8* 

100000.00 

15 

*7.69 

1*.19 

-2*.** 

.88 

16 

*7.77 

14.11 

-19.10 

.6* 

17 

*8.18 

1  *.03 

-17.19 

100000.00 
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5.  Axis  plot 

Axes  are  plotted  automatically  by  GSPP  in  connection  with 
profile  and  contour  plots,  unless  its  plotting  is  suppressed.  The 
subroutine  ACHSE  is  responsible  for  the  axis  plot.  It  is  designed 
for  virtually  all  different  cases:  arbitrary  direction,  scaling,  tick 
marks  right-or  leftbound,  scale  numbers  right  or  left  of  axis  in  four 
different  directions  (integer  multiples  of  90°),  variable  number  of 
decimal  digits,  variable  height  of  scale  numbers,  variable  tick  mark 
length,  variable  distance  (axis,  scale  numbers),  and  many  other 
options  more.  For  a  detailed  reference  see  the  comment  of  the  program 
listing. 

In  the  following  three  examples  of  axis  plots  are  shown.  Axis 
n  corresponds  to  the  n'  th  call  of  ACHSE  in  the  subsequently  listed 
test  program. 


PIGS  X5  ESsrc  QUMiiVY  ] 


FAXMPi  LS  OF  im'FEUENT  AMIS  PLOTS 

A  DESCRIPTION  OF  THE  INPUT  (AND  OUTPUT)  PARAMETERS  CAN  DE  FOUND 
IN  I  HI  SUBROUTINE  ' ACHSE  * 

LAI  L  f‘t  OTS  (O .  40.  -  i  w ) 

Pi  OT  UNIT  ; 


CALL.  A*: 
:£5  .UXv>.  i 
CALL  Ai 


TITSR  ()••!.  .  9,  -15.  ,  20 .  ,  5 .  .  2,0.  ,  5 .  ,  0 .  2 . 0 . 5 , 1 .  ,  1 .  ,  0 . 3 ,  i .  ,-1,0,1 
‘ITU) 

T!SF  ( J  >  •>.  .  1 1  .  190.  .  200 .  .  25.  ,2. AO.  .20.  ,0.25.  -0 . 4 , 1  .  ,  1  .  ,  0 . 25 . 1 


::  1.5.5 , 0  ,  n:<  ;) ,  NTH> 

:.M  L.  A  CHS  II  ( :■*  .  1  s .  .5  500 .  .  3500 .  .  500 .  ,  2 , 300 .  .  250 .  ,  ■  0 . 3  4  0 . 5 , 1  .  ,  1  .  ,  0 . 3 
£  5  >  ,  ••  1  ,  t  ,  J  ,  1  .  DX0J4TM) 

CALL  Pi  OT (h.  .0. ,999) 

S  I  OP 

t:  KM) 


% 


axis  No.  l 


Figure  5.1:  Examples  of  axis  plots 
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6.  Grid  plot 

Superimposing  a  grid  over  a  contour  plot  simplifies  a  possible 
interpolation  process.  GSPP  assures  that  the  individual  grid  interval 
is  identical  with  the  tick  mark  interval,  if  an  integrated  grid  plot  is 
requested.  But  a  grid  of  arbitrary  shape,  size,  and  orientation  can  be 
superimposed  also  externally.  This  is  done  by  the  subroutine  GRIPLO. 
A  number  of  options  can  be  chosen  by  the  user:  A  rectangular  grid 
can  be  full  line  grid,  a  dashed  line  grid,  or  simply  consist  of  open 
crosses  at  the  grid  points  (intersections  of  grid  lines);  the  grid  can  be 
arbitrarily  oriented,  can  be  shifted  arbitrarily  relative  to  the  lower  left 
comer  of  the  grid  area,  the  grid  distances  in  x-  and  y-  direction  are 
arbitrary,  the  plot  of  a  bounding  rectangle  (representing  the  grid  area) 
can  be  suppressed,  and  other  options.  For  a  detailed  information 

see  the  comment  statements  in  the  sub-routine  GRIPLO.  Some  grid 
examples  are  shown  below. 


I-XAMCUS  or  IHITTI.TNT  SRID  PLOTS 

;\  )>n;*.:UlPTI<lU  OF  INPUT  CARAME'l  ERS  CAN  DE  FOUND  IN  THE  SUBROUTINE 

'OR  1C!  O' 

CALL  Cl  OTSO> .  .0.  i  1«:>) 

CLOT  UNIT  :  CENTIMETER 


CAL  L 

ORICLOC  i 

•  1 

9 .  .  8 .  ,  6 .  .  (■» .  8 * 0 . 6 

,2. 

M 

,  4  ,  Q .  5 , 0 . 

,0. 

n 

CALL 

OKI  CLOU 

1 , 

,10.  j  8 1  ^  6 1  ^  0 »  ,  0 . 

|2*i 

4 ,  a .  2 ,  l . . 

•20 

CAI.  1 

CRIPLOC  1 

1  « 

•  ,1 

,7, 

0 .  1  , 1 . 5 , 

•  vl 

,0) 

CALL 

ORICLCK  i 

1  ] 

,,  1  .  .  8 .  ,  6 .  ,0.5)0. 

,2, 

2.5. 

3,  -0.2,2. 

,e. 

i) 

CALL 

■  o  r  co. . 

w, 

,999) 

O  f  tip 

. . . 
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7.  Title  plot 

The  description  of  plots  (bead,  axis  labels)  is  performed  by  the 
subroutine  TITLE,  a  part  of  GSPP.  The  main  purposes  of  this  program 
are: 

a)  find  the  number  of  lines  of  alphanumeric  characters  which  belong 
to  the  title  (in  connection  with  GSPP,  a  maximum  of  10  lines  is 
allowed;  in  isolated  applications,  a  maximum  of  100  lines  can  be 
plotted  by  a  single  call); 

b)  find  the  number  of  alphanumeric  characters  per  line  (max.  =80) 
and  the  maximum  number  of  characters; 

c)  find  the  maximum  actual  title  length  and  height,  compare  it  with 

the  corresponding  allowed  maxima,  and  reduce  the  symbol  height,  if 

(*  ) 

necessary. 

d)  Shift  the  title  lines;  four  different  title  line  patterns  can  be 
obtained:  the  title  lines  can  be  plotted  as  appearing  on  the  punched 
cards  (or  any  other  input  device);  the  title  lines  can  be  shifted  to 
the  left  (leftbound);  the  title  lines  can  be  shifted  to  the  right  (right- 
bound);  the  title  lines  can  be  centered. 

Moreover,  the  title  can  be  put  into  a  rectangular  frame,  called  title 
boundary  rectangle;  the  title  can  be  plotted  in  any  direction  (0°-360°), 

0°  is  the  horizontal  mode;  when  plotting  on  an  electrostatic  plotter 
like  the  Versatec,  the  line  width  can  vary  between  single  and  5-fold 
linewidth.  The  title  itself  is  stored  on  a  2-dimensional  array,  each 
column  representing  a  80-character  string  (1  card). 

Some  typical  examples  are  shown  below.  For  more  information 
(title  and  parameter  transfer)  see  the  comments  in  the  subroutine  TITLE. 


(*)  If  a  symbol  height  reduction  is  necessary,  a  message  will  be 
edited  on  the  assigned  output  unit. 


LXAHP*  fcS  OF  !'))  FFEUirNT  TITLE  PLOTS 
:! )  MENS  1  ON  1  )  T  C-JO,  10) 

COMMON  /111 1  I  PA/ Y ,  Si INK ,  TISHAX ,  XLTH .  ROT  ,  IPEN ,  I C ,  1 143 

r\  1TAUX/UII.  CMINCH,TH,T1TL.TH,  X  ,  ICY ,  1  0 ,MRL,  WAX  /IOSPAR/IO 

I  0! -6 

KLA!)  T 1  Hr.  NUMBER  Oi  TITLE  LINES 
RLVOKO.i:)  Nil 

Li  A!)  Hit  TIT  IF  PARAMETERS 

fvL.nl>  U. .  X)  X  ,  V  .  Sl-INR  .  THMA  X , UH , XL  TH ,  ROT  ,  CM  1  NCI  1 , 1 CV ,  I 143 , 1  PEN ,  I C ,  LC 
READ  TIT  DIM  (EACH  COLUMN  OF  TIT(...)  CORRESPONDS  TO  ONE  (1) 

irru  i  iNu 

UI.Al)(i-.S)  (Cirm  J) ,1-1 ,20) , J-1,NTL> 

:  Oiv’WA  i  (''OtVU 

CAI  t.  I  I  O  I  S  ('•> .  ,  0  .  .  1  H) 

CA1  I.  T  1  TL.K(  I  11,10) 

)  143-  <•> 
i  pi  n=  •;>. 

i\0  I .’*•)  i 

Y-'Y . 

CALL  I  1TI.ECI  IT,  10) 

1  C«1 
IPI.N'V 
ROT  =  w . 
y.  U. 

CALL  l  1TLEO  IT,  10) 

)  C ! '  0 
IT!  N~  1 
Y-o .  0 

CALL  i 1TLECT IT, 10) v 
CALL  i  I  O.T  ( 0 .  ,0.  , 9V9) 

OP 
I  NT) 


input  parameters: 

5 

0.,  20. ,  0.  3,  5.,  1.,  18.,  0.,  0.  3937,  0,  1,  5,  11,  0 
the  5  title  lines 


i>HT  uw»  u 
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TH1S  15  AN  EXAMPLE  OF  AM  AUTOMATICAL  TITLE  PLOT 
USING  THE  SUBROUTINE  TITLE  (A  PART  OF  GSPP) . 
THE  TITLE  CAN  HAVE  A  MAXIMUM  NUMBER  OF  100  LINES, 
IN  CONNECTION  WITH  A  PROFILE,  CONTOUR  OR  3-D  PLOT 

ONLY  10  LINES, 


f  *  ^  A  0 

"'S.4vS>  *"  »*' 

C0^ 


THIS  IS  AN  EXAMPLE  OF  AN  AUTOMATICAL  TITLE  PLOT 
USING  THE  SUBROUTINE  TITLE  (A  PART  OF  GSPP)  . 
THE  TITLE  CAN  HAVE  A  MAXIMUM  NUMBER  OF  100  LINES, 
IN  CONNECTION  WITH  A  PROFILE,  CONTOUR  OR  3-D  PLOT 

ONLY  10  LINES. 


THIS  IS  AN  EXAMPLE  OF  AN  AUTOMATICAL  TITLE  PLOT 
USING  THE  SUBROUTINE  TITLE  (A  PART  OF  GSPP). 

THE  TITLE  CAM  HAVE  A  MA.  ’ MUM  NUMBER  OF  100  LINES, 
IN  CONNECTION  WITH  A  PROFILE,  CONTOUR  OR  3-D  PLOT 
ONLY  10  LINES. _ _ 


Figure  7.1:  Examples  of  title  plots 

(The  title  and  the  top  correspond  to  the  first  call,  the 
title  at  the  bottom  to  the  last, call  of  the  subroutine  TITLE.) 
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8.  Range  division 

A  particularly  useful  element  out  of  GSPP,  is  the  subroutine 
RNGDIV;  it  is  frequently  used  in  GSPP  for  the  purpose  of  replacing 
the  human  decision  process  like  finding  "optimal"  contour  intervals, 
"optimal"  grid  distances,  or  "optimal"  tick  mark  distances.  "Optimal" 
is  interpreted  in  terms  of  a  most  reasonable  unbiased  decision. 

The  main  objectives  of  RNGDIV  are  as  follows:  given  a  range 
on  the  real  number  line  by  a  start-  and  endpoint  (lower  and  upper 
bound);  the  range  should  be  divided  into  intervals  of  the  form 

(0.1,  0.25,  0. 5)  *  10k 

with  k  an  integer;  the  maximum  number  of  intervals  can  be  chosen  by 
the  user;  there  is  the  possibility  of  adjusting  the  interval  start-  and 
endpoint  such  that  the  interval  start-  and  endpoints  have  "coordinates" 
of  the  form 

(a*  o.l,  B  *  0.25,  y*  0.5)  *  10* 

with  of  =  1,  . . ,  ,  ;  8  =  1,  2,  3  ;  y  -  l,  2.  The  adjustment  can  be  a 
range  extension  or  a  range  contraction  (e.g.  in  contouring).  The  program 
returns,  apart  from  the  calculated  interval  length  and  the  number  of 
intervals,  also  the  recommended  number  of  decimal  places  for  a  graphi¬ 
cal  representation  of  the  scale  numbers.  (Recommended  is:  3  significant 
digits,  if  the  maximum  scale  number  is  greater  than  100,  only  integer 
representation  of  the  real  scale  number.  )  A  special  application  is  the 
estimation  of  the  number  of  significant  digits  for  a  specified  range. 

In  the  sequal  a  couple  of  range  division  examples,  calculated  using 
RNGDIV,  are  listed.  Further  information  can  be  found  in  the  comment 
statements  to  RNGDIV. 


1'rUV  i-.-.G:- 

*  v  W  *  t 

L 


H 

c. 


IS  BI SCI ZUAJjin  FRACnCJBIS 
."^.ISffiSD  K>  r.'G  “ 

[£^f  -u.|  |ti£  jijo.w  j  T  Ii'jl  *  K.'JtiLi  i  '  •  »-'0.\  DcTaiLS  but  THH. 

tUt’ir  l  v  I  b  10  1 »  l  A  L  l  i  V 

x  npu i —0 
i‘.Ou  f-o 

viR  i  I  C.  I  miU  1  t  3  i 

kEAGi  iur-j  I  »"»  tti'S1*:  1  a«  *  Xu  *  NK  in  »  Ii*. 

CAlL  X.iouIV  l  Art  *Xt  »XL«  »  aLL  »L)AL»urtl.‘if  I  sM  »  R I  fOD  > 

Khl.,1  i'Jt'uT  uulHcT 

wt<  I  T  c.  (  ,1.J-J  l  »4  ,  A  A  ♦  ,s£  t  XL*  »  xLt  t  tiXL»  *010  r  I  l\l  n’jH  1  *.\.U 
uOl'O  1 

r  *jRii«  I'  l  l.i  »4A»»rLsT  Oh  SubROoTl.'IC.  RUoDI  V  ‘  r  dX» 

*'  V  /*.  At 

*c  t"  Xl  m\lrl  I  0  •'••'I 

r  0 1\  r.  h  1  l.i  r,  ♦  a  >  b  F  i  >.2  *4  i5  ) 

3  I  OP 
L  I*  Li 


XL A  XLL 

M  u  *  »  X  /  ) 


*  T  m«  u  u  x  V 

Or  bucKuC  I  1  iv r.  h ii 0 -  i  V 


X  A 


•ll«!  .b>' 
•i  !  2.  b(* 
-112.  3*3 


XE 


27 3. a  j 

2  7  3 . 4 '  ■ 
27  3.-*  • 


Xlh 


■  1  »  „  .  JO 

•  1 2  b .  v.  •: 
■12J.-J.1 


XL£ 


2  b*  • 

2  7  3  • 1 '  > 

4  7  3.0*1 


Dal 

NR  ill 

Iri 

OR  I 

OU 

23.00 

2'i 

-  I 

14 

-  1 

2  3.00 

21. 

1 

16 

-1 

2  3.0*"' 

0 

16 

-  1 

9.  Clipped  boundary  line  plot 

In  connection  with  contour  plots  in  regions  of  arbitrary  shape 
it  is  often  requested  to  plot  the  boundaries  of  these  regions  with  the 
restriction  that  the  boundary  lines  should  be  clipped  off  at  its  intersec¬ 
tions  with  the  boundary  of  the  rectangular  plotting  domain. 

The  subroutine  BNDPLO  is  designed  for  this  purpose.  It  accepts, 
in  principle,  an  arbitrary  number  of  (district)  boundary  lines  which  can 
be  plotted  with  different  line  widths  on  an  electrostatic  plotter  or  in 
different  colors  on  a  multi-color  plotter. 

In  order  to  clip  off  the  line  plot  at  the  boundary  of  the  rect¬ 
angular  domain,  each  line  element  passes  a  procedure  like  that  described 
in  section  2.8.5  which  provides  sufficient  information  about 
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the  location  of  the  line  element  relative  to  the  rectangular  domain.  The 
actual  intersection  uses  a  simple  line  intersection  algorithm.  If  so 
desired,  the  boundaries  can  be  mapped  according  to  an  arbitrary  mapping 
equation,  .x  =  x  (x,  y),  y  =  y  (x,  y),  which  has  to  be  provided  by  the  user 
(subroutine  GSPROJ). 

The  following  example  may  illustrate  the  use  of  the  subprogram. 

C  TEST  OF  THE  SUBROUTINE  'BNDPLO'.  FOR  DETAILS  SEE  THE 
C  COMMENTS  TO  'BNDPLO'. 

DIMENSION  F0RMD(19) 

COMMON  /XBOUND/X ( 1  99)  /YB0UND/Y(199)  /BOUI NF/NAE (29) 

*  /BOUPEN/NPEN< 19)  /CCLIP/XL , XU , YL , YU 

1NPUT=5 
NOUT-6 

C  READ  AND  ECHO  THE  NUMBER  OF  BOUNDARIES,  ITS  START-  AND 
C  EMDPOINT  LOCATIONS  ON  THE  VECTORS  X(.)  AND  Y(.)  AND  THE 
C  CORRESPONDING  PENUIDTHS  NPEN( . ) 

READ ( INPUT ,*)  NAEd) 

N=NAE( 1 ) 

NB=2*NAE  ( 1  )■*•  1 

READ! INPUT,*)  (NAE ( I) , 1-2, NB> 

READ (INPUT,*)  (NPEN( I ) , I"1 , N) 

WRITE (NOUT , 19) 

19  FORMAT (1H9,4X, 'BOUNDARY  LINE  INFORMATION  VECTOR  NAE (  .  )  '  , 

*/ ,5X , '  I  NAE (1 )  NAE(I+I)  NPEN( I /E) ' , //) 

DO  1  1*2, NB, E 
11=1+1 
12*1/2 

1  WRITE  (NOUT  ,11)  I ,  NAE(  I )  ,'NAE(  1 1 )  ,  NPEN  ( 12) 

N=NAE  (II) 

11  FORHATOH  ,4X, 15, 19, 18,119) 

C  READ  INPUT  FORMAT  FOR  BOUNDARY  COORDINATES 

READ( INPUT, 16)  FORMD 
16  FORMAT (19A4) 

C  READ  AND  ECHO  THE  BOUNDARY  COORDINATES 

WRITE (NOUT , 12) 

12  FORMAT (1H9,4X, 'BND#  I  X(I)  Yd)',//) 

N=NAE ( 1 ) 

DO  2  1=1, N 
1 S=NAE (2*1) 

I E=NAE (2*1 + 1 ) 

DO  2  J=IS, IE 

READdNPUT, FORMD)  X(J),Y(J) 

2  UR1 TE(NOUT , 13)  I  , J  ,  X ( J) , Y ( J) 

13  FORMAT ( 1H  , A X , 21 5 , 5X , 2F 1 9 . 2) 

C  READ  AND  ECHO  LOWER  AND  UPPER  X-  AND  Y-COORDINATES  OF 
C  THE  RECTANGULAR  PLOTTING  DOMAIN 

READdNPUT,*)  XL, XU, YL, YU 
UR1 TE (NOUT , 1  A)  XL, XU, YL, YU 

14  FORMAT (1H9,4X,'PL0TTIN6  DOMAIN  LIMITS  «',//, 5X, 

* ' XLOUER  ...  ',F19.2,5X,' XUPPER  ...  ' ,F19 . 2 , / , 3X , 

*'YLOUER  ...  ',F19.2,5X,'YUPPER  ...  ',F19.2) 

C  READ  AND  ECHO  PLOT  PARAMETERS  (ORISIN,  SCALE) 

READdNPUT,*)  X9,Y9,FAC 
URITE(NOUT, 15)  X9, Y9,FAC 

15  FORMAT (1H9,4X , 'PLOT  PARAMETERS  t',//,3X,'X9  ...  , 

*FI9.2,5X,'Y9  ...  '  ,F19.2,5X, 'FAC  ...  ',F19.2> 

C  INITIALIZE  PLOT 

CALL  PLOTS (9 ,  ,9.  ,  19) 

CALL  FRAHE((*  ,9.  ,199.  ,9.  ,79.  ) 

C  PLOT  THE  BOUNDARIES 

CALL  BNDPLO (X0,Y9,FAC,9,9,9) 

C  PLOT  THE  BOUNDARY  RECTAN6LE  OF  THE  PLOTTING  DOMAIN 

CALL  RECT(Yw, XO,  (YU-YD/FAC,  (XU- XL) /FAC, 9.  ,  1 ) 

C  STOP  PLOT 

CALL  PLOT (9. ,9. ,999) 

STOP 

END 
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Input  data  to  the  above  program 

15 


e 

,19,86,87, 35, 

36,49,31 , 

3 

1,2 

,3,1, 

2 

4 

m 

,EF10 

.2) 

3 

i 

E 

45.89 

-El. 89 

6 

l 

3 

41.69 

-16.49 

7 

l 

4 

41.89 

-9.49 

• 

l 

5 

45.89 

-15.29 

9 

l 

6 

59.69 

-14.69 

1* 

i 

7 

54.99 

-17.29 

11 

i 

8 

59.29 

-18.99 

IE 

l 

9 

59.99 

-18.89 

13 

l 

19 

45.89 

-El. 89 

14 

E 

14 

55.99 

5.E9 

15 

2 

15 

62.99 

1.89 

16 

2 

16 

68.99 

-4.69 

17 

2 

17 

84.69 

-6.99 

18 

2 

18 

86.49 

-3.89 

19 

2 

19 

88.29 

-9.89 

Eft 

2 

Eft 

88.29 

1.49 

El 

2 

21 

81.89 

6.49 

EE 

2 

22 

78.99 

11.69 

83 

E 

E3 

78.99 

El  .Eft 

24 

2 

24 

75.69 

E4.99 

E5 

2 

25 

62.49 

E2 . 6ft 

E6 

8 

E6 

50.49 

E5.E9 

27 

2 

27 

45.89 

23.69 

E8 

E 

E8 

46.99 

12.89 

E9 

E 

E9 

49.89 

7.80 

3ft 

2 

3ft 

55.90 

5.29 

31 

3 

36 

71.49 

11.89 

38 

3 

37 

75.20 

4.89 

33 

3 

38 

62.90 

6.40 

34 

3 

39 

55 . 29 

9.89 

35 

3 

4ft 

54 . 90 

12.  *9 

36 

3 

41 

54 . 90 

16.4ft 

37 

3 

42 

64.80 

16.2ft 

38 

3 

43 

67.20 

19.29 

39 

3 

44 

71.40 

11.89 

4ft 

4 

45 

39.49 

35.69 

41 

4 

46 

54.80 

32 . 99 

42 

4 

47 

55.90 

49.29 

43 

4 

48 

50.60 

35,49 

44 

4 

49 

39.40 

35.6ft 

45 

5 

52 

193.49 

22.69 

46 

5 

53 

84.80 

33.49 

47 

5 

54 

68.99 

28.99 

48 

5 

55 

61.80 

29.89 

49 

5 

56 

63.99 

35.99 

5ft 

5 

57 

79.69 

39.8ft 

51 

5 

58 

193.49 

22. 6ft 

52 

25. 

E,  92 , 

2, -12.5 

40.9 

53 

3., 

3* 

Figure  9.1  Example  of  clipped  boundary  line  plot  corresponding  to 
the  above  input  data 
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